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ABSTRACT 

We map the distribution of dust in M31 at 25 pc resolution, using stellar photometry from the 
Panchromatic Hubble Andromeda Treasury survey. The map is derived with a new technique that 
models the near-infrared color-magnitude diagram (CMD) of red giant branch (RGB) stars. The 
model CMDs combine an unreddened foreground of RGB stars with a reddened background population 
viewed through a log-normal column density distribution of dust. Pits to the model constrain the 
median extinction, the width of the extinction distribution, and the fraction of reddened stars in 
each 25 pc cell. The resulting extinction map has a factor of >4 times better resolution than maps 
of dust emission, while providing a more direct measurement of the dust column. There is superb 
mo rphological agreemen t between the new map and maps of the extinctio n inferred from dust emission 


by Draine et al. (2014). However, the widely-used Draine & Li (2007) dust models overpredict the 


grains are significantly more emissive than assumed in Draine et al. ( 20I4|). The observed factor of ^2.5 

discrepancy is consistent with similar findings in the Milky Way by 

Flank Collaboration et al. (2014), 

but we find a more complex dependence on parameters from the 

Draine Li (2007|) dust models. 

We also show that the discrepancy with the Draine et al. (2014) 

map is lowest where the current 


interstellar radiation field has a harder spectrum than average, we discuss possible improvements to 
the CMD dust mapping technique, and explore further applications in both M31 and other galaxies. 
Subject headings: ISM: dust, extinction, ISM: structure, galaxies: ISM, galaxies: stellar content, 
galaxies: structure 
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1. INTRODUCTION 

Dust plays an increasingly important role in extra- 
galactic astronomy. It has long been known that dust 
shapes the observational properties of disk galaxies, par¬ 
ticularly in the optical and ultraviolet (e.g., Disney et al. 
19891 IXilouris et al.||19991 |Misiriotis et al. 200 1[ |Fierini 


et al 


.]2QQ4r. 

Driver et al. |2007| |Bianchi||20b8| jPopescu et al.|[2011)^ 

However, over recent decades, dust has become a wide 


Tuffs et al 


20041 IMollenhoff et al j(2QQ6 


spread source of st udy in its own right (Savage & Mathis 
1979 Draine|2003), thanks in large part to observational 
facilities that directly probe emission from the dust in 
the mid- and far-infrared, and at sub-millimeter wave¬ 
lengths. This emission has also become a widely used 
tracer of key astrophysical quanitites, including the star 
formation rate and the interstellar radiation field (see 
review by Kennicutt & Evans||2012 ). 

Dust has also emerged as one of the most effective trac¬ 
ers of cold, dense gas. Molecular gas and cold Hi are the 
immediate precursors to star formation (see, for example 
Bergin fc Tafalla1|2Q071 |McKee fc Ostriker||2QQ7[ |Hen- 


nebelle & Falgarone|20i2l), and the properties of this gas 


is likely to be coupled to the ability of the interstellar 
medium (ISM) to form stars. Unfortunately, these cold 
gas components are extremely difficult to trace. Cold 
molecular hydrogen is nearly impossible to see in emis¬ 
sion, and cold atomic Hi cannot be reliably distinguished 
from warmer phases without absorption line techniques. 

The study of cold, dense gas has been fundamen¬ 
tally change d by the widespread use of dust extinction 


as a tracer (Lilley 1955 Dickman 1978 Frerking et al 
1982). Methods tnat used star counts to identit y dust 
extinction had been in use for many years (e.g. , |Dick- 


1997 


.V_/OXV_/XJ. XXJ. J.V_/X ^ 

19781 ICernicharo fc Bachiller|1984 Cambr4sy et al 


Arce &: Goodman||1999[ [Dobashi 


et al.| ^Q5[ among many others), but were gradually 


supplanted by new methods taking advantage of the 
growing availability of near-infr ared imaging, particu- 


larly due to the all-sky 2MASS ( Skrutskie et al. 2006) 
and the UKIDSS/Galactic Plane (Lucas et al.|2UU8|) sur- 
veys. Many groups have developed optimized method- 


magnitude) diagrams (e.g., Lada et al.|1994| Giardi et al. 

1998 

ILombardi & Alves 2001| [Cambresy et al.[ 2002; 

Lorn 

)ardi| 20051 IFroebrich & del Burgo 2006| Lombardi 

2009 

Gonzalez et al.|201ip to identify molecular clouds in 

the 1\ 

dilky Way and to map their structure ([Giardi et al. 

1998 

Alves et al.|1998t[Lada et al.|19991 [Alves et al.|2001[ 

Gam 

3resy et al. 20021 Teixeira et al.|2005 Lombardi et al. 

2006 

Froebrich et al.||2007| ILombardi et al. 2008| Kain- 

ulainen et al.||2009| IRowles & Froebrich[|2009 Lombardi 

et al.|2010[ IRoman-^uhiga et al.||2010| IFrieswijk & Ship- 

man||2010| ^candariato et al.|2011 Dobashi||2011t Kain- 

ulainen et al.||2011a 

bl Gonzalez et al. [[20121 [Alves et al. 

2014| Schlahy et al. 

2015| and many others), with even 


even 

denser clouds (Vasyunina et al.[ 

[20091 

[Butler & Tan 

2009[ 

i IMajewski et al.||20111 IButler & Tan| 

POTW 

|JN idever 

et al, 

.[20121 Kainulainen & Tan 2013), and otl 

ler work 


et al. 2009). The resulting maps are likely to be more 
direct tracers of the total column density t han meth- 


ods based on dust or molecular emission (e.g., Goodman 


et al.||2009). 

The impact of extinction mapping can easily be seen 
in the diversity of problems they have been used to ad¬ 
dress. These maps have been used: to measur e the 
statistics of the column density distribution (e.g 


et al.|19^|Kainulainen et al.|2009l Lombardi et al.|^i^010[ 


||2014D __ _ . 

et al.B20091 IRowles fc Froebrich||2011[ [Kainulainen et al.' 

2011b[ ~ ?tutz &: Kainulainen|2015D; to explore the origins 


Froebrich fc Rowles||2010| [Schneider et al.||201H [Alves 
et al. 1120141) and its connection to star forrnation dLada 


of 'Taiion' s Laws” (|Kauflniann et al. 2010 Lombardi 
et al. 112010 Beaumont et al.f 2012| Ballesteros-Paredes 
et al.|2012 ); to meas ure the cloud structure function and 


clump statistics (e.g., Padoan et al.[2002[ 200^ Kirk et al 


2006 Lombardi et al. 2008 ' 2010 ' Kauttmann et al.|l 
to stimy individual molecular cores and globules (e.g 
Teixe ira et al.||2005l IRoman-Ziihiga et al.||20091 |Racc r 
et al 



molecular gas (e.g., Lada et al. [1 

994 

: Hayakawa et al. 

20011 Harjunpaa et al.|2004||LomD 

ardi et al. 120061 [Kain- 

ulainen et al. 2006 Pineda et al.||2010|); to distinguish 

among models of turbulence ([Padoan et al.[[ 19971 

[Feder-[ 


rath et al.[[2010 

Kainulainen i 

T'an[[2013 

l; to derive the 

3-dimensional structure of the ISM (e.g.. 

Sale & Magor- 

rian|2014 Schlafly et al. 12015 

Green et a 

.|2U15|); and to 

constrain models of the dust itself (e.g., Roy et al.|2013). 


the study of external galaxies. The Magellanic Clouds 
are sufficiently close to apply similar map ping techniques 
using NIR colors of individual stars (e.g., [Dobashi etaT 


2008 2009), although their low extinctions a llows optical 


and ultraviolet data to be used as well (e.g., Harris et al. 


1997 Zaritsky et al. 2002 2004). At larger distances 


where individual stars can no longer be resolved, extinc 
tion maps are typically based on modeling the spectral 


energ y distribution in individual pixels (e.g., Zibetti et al. 


2009), but usually only as a way to removing the ettects 


of dust, rather than as a route to studying the cold ISM 
itself. 

In this paper, we focus on using M31 to bridge these 
two regimes. Andromeda is close enough that we can 
probe the ISM on the scales of molecular clouds. How¬ 
ever, thanks to the large area covered by th e Panchro¬ 
matic Hubb le Andromeda Treasury (PHAT; jPalcantcm 


et al. 2012), we can also cover large areas, allowing us 


to view the entire cold ISM, rather than just individ¬ 
ual clouds. This approach gives us our first view of the 
statistical properties of molecular clouds across a large 
fraction of a massive spiral galaxy. 

Thanks to its proximity, M31 has been the target of 
many previous dust studies, including surveys of dust 
emission using space-bas ed mid- and far-infrar ed (FIR) 
imagers on board IRAS ( Deyere ux et alj|1994 ), Spitzer 


(Barmby et alJ| 2Q^ Gordon^^. 2006p , and Herschel 


(4ritz et al. [Groves et al. T'hese observations 


gas ratio, and dust composition (e.g.. Ford et al.[ 

2013 

Montalto et al. [20091 Tabatabaei & Berkhuijsen| 

OT 

Leroy et al.|20111 

Groves et al.[[2012| [Smith et al.[ 

2012 

Draine et al.[[2014 

), particularly when complemented by 


1986| [Braun et aL][2QQ9| [Ghemin et al.|[2UU9[) and GO 
dblieten et al.[[2UU6| Rosolowsky et al.|[2U07[ Tbsaki et al. 
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2007 though the latter two cover limited area). 

Some of these studies have used the dust emission and 


other trace r s to derive extincti on maps (e.g., [Montalto 


et al 


20091 [Praine et al. 2014). Most recently, [Draine 


et al. HMr used mid- and tar-infrared data to derive 
dust column densities throughout M31. These maps use 
state of the art models of dust to derive the spatial dis¬ 
tribution of dust composition, heating, and dust column 
density. The models also make specific predictions for 
the extinction within M31. 

In this paper we introduce a technique for probing the 
cold dusty ISM. We take advantage of t he Panchromatic 


Hubb le Andromeda Treasury’s (THAT; Dalcanton et al. 


2012) large database of NIR HST observations, and use 


photometry of individual RGB stars to derive the distri¬ 
bution of dust extinction on 25 pc scales. Specifically, 
we use the structure of the RGB seen in NIR color- 
magnitude diagrams (CMDs) to fit for the distribution 
of extinctions along the line of sight. The methodology 
therefore gives us both the median extinction and the 
width of the extinction distribution in each resolution 
element (oixel). 

In Sec. we give an overview of this technique, and ex¬ 
plain its connection to observations of molecular clouds 
in the Milky Way. In Sec.j^we derive the expected distri¬ 
bution of color and/or reddening for a log-normal distri¬ 
bution of dust embedded in thicker stellar disk. In Sec. 
1^ we give details of how we fit the parameters of the 
dust+star model to data from the PH AT survey, and 
discuss the distribution and accuracy of the derived pa¬ 
rameters in Sec. We show the global dust extinction 
map and compare it to the extinction inferred from fits 
to dust emission Sec. (H We then conclude in Sec. 0 


2. OVERVIEW OF MEASURING EXTINCTION WITH CMDS 

Mapping extinction in the Milky Way requires cop¬ 
ing with the complex distribution of both dust and stars 
along the line of sight. In external galaxies, however, the 
spread in distance is negligible, and all stars can be as¬ 
sumed to be at the same distance. This difference allows 
the full color and magnitude distribution of stars to be 
used when constructing extinction maps. 

The impact of dust can be most easily perceived in the 
NIR GMD, which is populated almost entirely with RGB 
stars that occupy a small range of intrinsic color at given 
magnitude. Theoretically, the RGB should form a nar¬ 
row, nearly vertical sequence in a NIR GMD. Empirically, 
however, the GMD structure of the RGB is complex and 
spatially varying, as shown in Figure In many regions 
of M31, the NIR RGB is quite thin (as seen later in Fig- 
ure|^, but in other regions it appears broad, and in some 
cases even bimodal. The broadening and/or bimodality 
is always seen redward of the thinner, bluer component, 
the latter of which always stays in essentially the same 
place on the GMD. Given that the color distribution of 
the bimodal RGB is not sensible in any stellar popula¬ 
tion model, and that the unusual morphologies change 
rapidly over very short distances, the natural explana¬ 
tion for the effects seen in Figure 0 is dust. 

The GMD morphology in Figure]l]can be explained as 
follows. The thinner, bluer RGB sequence is due to old 
or intermediate age stars seen in front of the cool dusty 
ISM. The redder stars are then those that lie behind the 


layer of dusfp^ The cool ISM is highly structured, and 
thus we see rapid small-scale spatial variations in the 
behavior of the reddened stars. Based upon this inter¬ 
pretation, we developed a method to simultaneously fit 
the unreddened and reddened stars to derive the distri¬ 
bution of dust extinction, which we describe in detail in 
Secs. p&Hl 

At its heart, the method is based on recognizing that 
the RGB stars are point-like “samples” of the column 
density distribution function of the dust. We demon¬ 
strate this effect in Figure The first panel shows 
a 2MASS-based extinction map of the Ori on molecular 


cloud, kindly provided by J. Kainulainen (Kainulainen 
et al. 2009^^ We have adjusted the angular scale of 
the niap to oe consistent with being located at the dis¬ 
tance of M31. The grid on the image indicates regions 
25 pc across, which is the scale we use in our analysis 
below. The second panel shows the distribution of ex¬ 
tinctions within each of these grid regions. In the PHAT 
NIR photometry, there are typically ^100 RGB stars 
that fall within a region this size, and 20-80% of them 
he behind the dust. The RGB stars therefore sample the 
dust structure on scales below that of the 25 pc grid, and 
even adjacent RGB stars may fall on regions with quite 
different extinctions. 

If we assume that any stars coincident with regions 
in Figure are primarily found in a stellar disk that is 
significantly thicker than the cool ISM (consistent with 
the h igh velocity dispersi on of ^90kms“^ seen in RGB 
stars; Dorman et al.|20I5), some fraction of the stars will 
be behind the dust layer, with the exact fraction depend¬ 
ing on the disk geometry. These stars will be reddened, 
and will sample the local distribution of reddening in the 
molecular cloud. By analyzing the reddening distribu¬ 
tion, we can infer the statistical properties of the extinc¬ 
tion distribution of the gas layer that they lie behind, for 
a given choice of extinction law. The same approach ap¬ 
plies for any dusty gas layer, not just the densest regions 
identified as molecular clouds, as long as the amount of 
reddening is large enough to produce a measurable shift 
in the RGB. 

As an example, the third panel of Figure shows the 
mean extinction calculated in 25 pc bins (i.e, the same 
width as the grid lines in the upper panels), if the same 
region of the Orion molecular cloud shown in the upper 
panels were observed in M31; this bin size matches the 
fiducial bin size used in our analysis below, and provides 
a good balance between ensuring a large number of stars 
per bin, and maximizing resolution. To improve the spa¬ 
tial sampling of the distribution, we calculate the mean 
in 4 different grids, each shifted by half of a 25 pc pixej^ 
and then interleave them into a grid with 12.5 pc pixeL 


Based on the Milky Way (e.g.,|Combes|l99l||Juric et al.|2008|), 
the cold ISM is almost certainly in a tJiinner layer tJian tJie RUB 
stars, which in MSI are primarily in a > 0.5 kpc thick disk (as we 
show in a companion paper). Thus we can safely neglect the effects 
of stars that are embedded within the dust. This assumption would 
not necessarily hold when analyzing a much younger population of 
stars that is expected to be largely embedded within the cool ISM. 

Although only a small fraction of MSI’s ISM is likely to be 
in the form of molecular clouds comparable to Orion, this example 
still serves as a useful illustration of the technique. 

Although not an exact equivalent, this approach is close to 
that of dithering to produced a Nyquist sampled image — i.e., we 
generate 12.5 pc pixels of a map with native 25 pc resolution. 
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Fig. 1.— The spatial variation in the NIR color magnitude diagrams, shown in 4 adjacent 20^' bins within the disk of M31. Each CMD 
is dominated by red giant branch stars. There is a narrow, bluer RGB sequence that is similar in all regions, but that is accompanied 
by a highly variable redder RGB sequence, which varies in both width and position. RGB stars trace a sufficiently old stellar population 
that they should be well mixed at these spatial scales. The spatial variation in the redder population is therefore best explained by spatial 
variation in the dust that is obscuring some fraction of the RGB stars. 




Fig. 2.— Simulation of sa mpling the extinction distr ibution of the Orion molecular cloud at the distance of M31. (First Panel) The 
extinction map of Orion from |Kainulainen et al.| (|2009|), with the angular scale adjusted to place it at the distance of M31, and the RA 
and Dec shifted to a location m lVi3Fs iUkpc star-torrning ring. Grid lines are 25pc (= 6.645'' apart). (Second Panel) The distribution 
of extinctions within each grid region (yellow), superimposed on the map from the top left. The same range of extinction is shown in 
each panel (0 < Ay < 5 mag). The distribution of extinctions can be reasonably characterized as a log-normal distribution in most cases. 
(Third Panel) The mean extinction, calculated using a dithered grid of 25 pc, in steps of 12.5 pc. In spite of the low resolution, the basic 
structure of the molecular cloud is still clearly visible. (Fourth Panel) The mean extinction in the same grid as in the lower left panel, 
but now calculated using only ^^25 randomly sampled pixels from the high resolution left-hand image in each 25 pc bin. Even this modest 
number of samples is sufficient to capture the same morphology seen in the middle panel, with a slight tendency to miss some of the highest 
extinctions. 


The resulting map of the mean extinction clearly traces 
the same structure as the actual molecular cloud, albeit 
with lower resolution. The dynamic range is reduced be¬ 
cause we are tracing the mean, and not the extremes of 
the distribution, and are convolving the data over the 
scale of a pixel. 

The fourth panel of Figure shows the same calcula¬ 
tion of the mean extinction, but instead of averaging all 
the points in the original image, we use only averages 
from randomly selected lines of sight, which were laid 
down with an average spatial density of ^ 25 samples 
per 25 pc bin (i.e., equivalent to what would be observed 


for an equivalent number of stars observed through the 
cloud). In spite of the fact that the fourth panel has dra¬ 
matically fewer samples per bin than in the third panel, 
the map looks nearly identical. Thus, even a relatively 
coarse sampling of the extinction distribution is sufficient 
to map out broad statistical properties of the reddened 
distribution. 

The example in Figure]^ bolsters our assumption that 
we can neglect depth effects within molecular clouds. 
The features in the Orion map have typical sizes that 
are less than 50 pc in their longest extents, and that 
are much smaller perpendicularly. Assuming that these 
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small projected sizes are also comparable to the depth 
of the cloud along the line of sight, then the cloud must 
have little depth compared to the overall distribution of 
RGB stars. The scale height of M31’s stars is of or¬ 
der ^ 0.5 — 1 kpc, which is more than a factor of 10 
times larger than the likely depth of the cloud. Thus, 
the distribution of extinctions for groups of neighboring 
stars primarily reflects sampling the distribution of ex¬ 
tinctions across the face of individual clouds, rather than 
a distribution of line-of-sight extinctions from RGB stars 
at different depths. It also suggests that when modeling 
the extinctions of comparable regions, we can most likely 
neglect any depth effects, and treat the dust as a thin 
screen within the stellar disk, such that the vast major¬ 
ity of stars are either in front or behind the dust layer. 
This assumption is likely to also hold even if there are 
multiple molecular clouds along the line of sight, given 
the typical thinness of molecular gas distributions. 

That said, the neglect of depth effects may be less valid 
in regions where most of the gas is in the atomic phase, 
rather than in the dense molecular phase that has been 
the typical domain of dust extinction studies in the Milky 
Way. M31 is Hi dominated, and thus in many regions 
the gas layer may well be intrinsically thicker than is 
characteristic for individual molecular clouds, given Hi’s 
typically larger velocity disperions. In practice, however, 
it seems unlikely that most of the Hi is in a layer as thick 
as the stars. The Hi dominated regions are also likely to 
be the lower extinction regions where our methodol ogy 


2012| and references therein) suggest that ISM densi- 
ties follow a log-normal distribution, both in terms of 
the volume density and the projected surface density. 
Observationally, anal yses of dense molecular clouds and 
the diffuse ISM (e.g., [Berkhuijsen fc Fletcher 2008 Hill 


et al. 


2Q Q8| have found that the probability distribution 


function of dust and/or gas follows a broad distribution 
around a peak, that in most cases appears to be well fit 
by a log-normal density distribution. Although there are 
recent questions about whether the distributions are best 
fit by a log-normal or by broken power laws (particularly 
at low de nsities where the measurements are most uncer - 
tain; e.g., Alves et aL]|2Q14 Lombardi et al.||2Q14[|2Q15 ) 


or whether deviations from log-normal distributions are 
the result of integrating over lar ge areas that sample very 
different cloud conditions (e.g., Brunt|[2M5 ), or whether 


somewhat different fitting functions are p referred in nu¬ 
merical simulations (e.g., [Hopkins 2013), these broad. 


peaked distributions appear to be a generic feature of 
the ISM, and are not tied to any particular gas phase. 

The only consistent deviations from a log-normal den¬ 
sity that have been noted are in extinction studies of 
the very densest regions, where some molecular clouds - 
particularly those that are star-forming - are observed to 


is most lik ely to break down for other reasons (Sec. 5.3 
and 5.3. 2| below). 


Before proceeding, it is worth noting the similarities 
and differences with other extinction mapping efforts. 
By concentrating on mapping a large area of an external 
galaxy, our work is similar in spirit to the series of Har 


(i.e., to high extinctions; Kainulainen et al.||2011b 

Rus- 

seil et al. 112013 

Schneider et al. 

PJT3 

Alves de O 

iiveira 

et al. 2014 He 

hneider et al. 20] 

L5| Al 

breu-Vicente et al. 


Tassis et al.|2QlQl|Ballesteros-Paredes et al.|2Qlll|Kritsuk| 


et al.|2QfT ). However, the fraction of a cloud’s rnass that 


is in this power law tail is typically small, such that the 
majority of the area and/or volume of the cloud samples 
the log-normal distribution. For example, in the map of 
Orion above, no more than 5% of the area falls in the 


et al.||1997| 

Zaritsky et al.||2002, 2004). However, we use 

the 25 pc analysis bins. 


only 2 ^IR 

hlters rather than fitting the full spectral en- 

Given its ubiquity, we 

expect the log-normal distribu- 


ergy distribution of the stars, and limit our analysis to 
only stars on the RGB. These choices reduce the data 
requirements and reliance on stellar models, but at the 
expense of being unable to probe differences in the dust 
column experienced by different kinds of stars. In addi¬ 
tion, our interpretation of the reddening distribu tion is 
sampling of a log-normal distribution”, whereas jZarit-j 


sky et al. (2004) largely interprets the width of the extinc- 
tion distribution as being due to stars’ different depths 
within a thicker diffuse d ust layer, although with some 
dumpiness being present (Harris et al.j 1997). Likewise, 
our approach has similarities to the IVIilky Way extinc¬ 
tion mapping described above, given that they are both 
based in modeling the NIR color distribution. However, 
we are able to use magnitude information as well as color 
information, thanks to the lack of any significant distance 
spread along the line of sight. This allows us to use only 
two filters (i.e., one color), whereas Milky Way studies 
typically employ 3 (i.e., the J, iL, and Kg filters of the 
2MASS survey). In an appendix, we discuss possible ex¬ 
tensions of this technique to bluer filters (Sec.|^, which 
could in principle be more sensitive to lower extinctions. 

3. EXPECTED DISTRIBUTION OF Ay 

Simulations of ISM turbulence (see reviews by 
Ballesteros-Paredes et al.pQllj [Hennebelle fc Falgaronej 


tion of ISM densities to be imprinted on the distribution 
of reddenings seen in the color-magnitude diagram. Since 
these reddened background stars are simply sampling dif¬ 
ferent lines of sight, they should experience the same dis¬ 
tribution of ISM column densities seen in the Milky Way, 
which will then be manifested as a log-normal distribu¬ 
tion of reddening. The one caveat is that the extinction 
distributions reported in the Milky Way are typically cal¬ 
culated over size scales that may not agree with our fidu¬ 
cial 25 pc bins. In cases where we analyze smaller areas 
we might expect less averaging over ISM properties, and 
thus potentially larger departures from the log-normal 
assumption. This effect can be seen in the second panel 
of Figure where three of the sub-regions show more 
complexity in their extinction distributions than can be 
easily described by a log-normal. Even in such cases, 
however, the log-normal still can provide a meaningful 
functional form for characterizing the typical extinction 
and its dispersion. 

Based on the above, for the remainder of our analysis 
we model the effect of dust reddening on the predicted 
distribution of Ay as follows. We first assume that all the 
extinction in an analysis region comes from a single layer 
of dusty gas along the line of sight, across which there 
is a log-normal distribution of line-of-sight values of Ay. 
Moreover, we assume that the gas layer has negligible 
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thickness compared to the stellar distribution, allowing 
us to assume that stars either experience all of the dust 
column, or none of it. The probability distribution pa 
for Ay then becomes a (5-function at Ay = 0 normalized 
to the fraction of stars that he in front of the gas, plus a 
log-normal function PA,gas describing the distribution of 
extinctions across the face of the gas layer, normalized 
to fred^ the fraction of stars that lie behind the layer: 


PA{Ay\fred, Ay , Cf) dAy = (1 - fred) ^{Ay) dAy 

+fred PA,gas{Ay\Ay,a)dAy, 

where the extinction distribution of the dusty gas layer 

Pa, gas is 


PA,gas{Ay\Ay,a) dAy 


(In (Ay / ^y)) 


Ay\/27ra‘^ 

The parameters of the log-normal distribution PA,gas 


dAy. 

( 2 ) 


Ay^ which is the median extinction of the stars behind 
the dust, and cr, which is a dimensionless parameter that 
sets the width and skewness of the log-normal distribu¬ 
tion; smaller values of a produce more symmetric distri¬ 
butions with less significant tails to large values. These 
quantities can be related to the mean extinction {Ay) 
and the standard deviation of the extinction a a using 
equations found below in Sec. |3.1[ In our model, the 
median extinction Ay is independent of the fraction of 
stars that are actually reddened, and will therefore have 
the same value whether the fraction of the stars behind 
a particular molecular cloud is 10% or 90%. This treat¬ 
ment could be easily generalized to multiple gas layers 
along the line of sight, but in practice we expect this 
occurrence to be rare. 

When modeling the CMD, the majority of the leverage 
on measuring pa does not come from measuring extinc¬ 
tion, but instead from measuring the reddening in the 
NIR FI low — F160W coloiF^ We therefore transform 
the extinction probability distribution into one of red¬ 
dening for an arbitrary color X — Y . Once the values 
of Ax/Ay and Ay/Ay are fixed by adopting a specific 
attenuation curve Ax/Ay, the reddening E{X — Y) be¬ 
comes 


E{X -Y)=Ay 


Ay 




( 3 ) 


Changing variables from Ay to the reddening £xy = 
E{X — Y), the probability distribution of reddenings be¬ 
comes 


PS,gas {£xY\A.v,<y,A\/^v)d£xY = 

£xy 


Pa, gas {Ay 


Ax 
Av ■ 


A 
Av 

■^y 


j\Ay,a) 


( 5 ) 


dA 


y. 


For the specific case of NIR data in the PHAT data 
set, we adopt Ap^^ow/Ay = 0.3266 and Apieow/Ay = 
0.2029 from [Girardi et al.| (|2008|)’s a pplica tion of a 
standard Ry = 3.1 Cardelli et al. ( 1989| ) attenua¬ 


tion law to typical RGB spectrum with a NIR color 
of FI low — F160W = 0.75. These ratios vary some¬ 
what with the temperature of the star, and thus are 
not constant for all of the s t ars on the RGB. However, 
fitting to the |Girardi et al.| (|2008|) isochrones for RGB 
stars shows that the color iepenc 


rat the color 
with A FiiQ\y / Ay = 0.3448 


ependence is quite weak, 
0.0243(F110W - F160W) 
and Aneow/Ay = 0.2061 - 0.0043(F110W - F160W). 
The RGB stars in our analysis have intrinsic colors in 
the range 0.5 < FllOW — F160W < 1.0, and thus 
the adopted extinctions are good to ±2% in FllOW 
and ±0.5% in F160W throughout the RGB, assuming 
Ry = 3.1. These scalings would be ^15% larger if the 
true attenuation law is much steeper {Ry = 5), with 
Apiiow/Ay increasing by 0.051, and Afiqqw /Ay by 
0.031. This change would be smaller than our typical 
precision. 

We will quote results in terms of Ay to simplify com¬ 
parisons with previous work and theoretical calculations, 
but the reader should note that the fit for Ay is driven 
almost entirely by F(F110W — F160W) = 0.124Ay. In 
addition, our adopted ratio F(F110W — F160W)/Ay 
is 6% smaller than the canonical value for a solar-type 
star, due to the underlying cooler temperature we have 
assumed. Although the exact conversion from reddening 
to extinction depends on assumptions for the underlying 
stellar spectrum and for the attenuation law, the former 
does not vary dramatically among RGB stars, and the 
latter appears to be quite stable in the NIR among dif- 
ferent lines of sight in the Milky Way (e.g., Martin & 
Whittetll^ . '-' 


3.1. Other Useful Quantities 

When interpreting dust maps, it is frequently useful to 
use the mean extinction (Ay), which is a better measure 
of the integrated column density of gas along the line 
of sight for a fixed dust-to-gas ratio. For a log-normal 
distribution, the mean depends on both the median Ay 
and the dimensionless width a: 


Ps{£xY\fred, Ay Ax/Ay) = (1 - fred) S{£xy) 

Yfred PS,gas{£xY\Ay, (7, Ax/Ay), 

where 

23 While extinction does offer some constraints through the 
structure of the RGB luminosity function, the intrinsic narrow¬ 
ness of the RGB makes color a far more sensitive indicator of the 
presence of dust. Other features — such as the tip of the RGB 
— have the needed intrinsic narrowness in luminosity, but do not 
contain many stars, making them a noisier probe of extinction. 


{Ay)=Ave^"/\ (6) 

In some circumstances it may be useful to consider the 
actual width of the log-normal in magnitudes of extinc¬ 
tion (i.e. the standard deviation cr^), rather than the 
dimensionless parameter a. If needed, one can relate a 
to the standard deviation a a of the distribution using 
either 

g\ = Ay e^ {e^ — 1 ) 


( 7 ) 
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or 


7^ = In 


1 + Vl + 4u;2 


( 8 ) 


where w = cjaI^v- Note that a a scales linearly with 
the median Ay^ such that in regions of high extinction, 
the same value of a corresponds to an intrinsically wider 
distribution. Because of this dependence, we choose to 
fit for the dimensionless width a, which is more likely to 
be invariant across the disk. 

4. MODELING THE OBSERVED DISTRIBUTION OF 
REDDENING 

The basic outline of our fitting for the distribution of 
reddening (i.e., ps) is as follows. We first use NIR data 
from the PHAT survey to generate high- qual ity CMDs, 
focused on the well-populated RGB (Sec. 4.1). We then 
create an empirical model of the unreddened foreground 
as a function of position, by isolating small subregions 
with extremely low reddening, as indicated by their nar¬ 
row, blue RGB and low FIR emission (Sec. |4.2| ). Assum¬ 
ing (1) that the RGB stars are well-mixed along the line 
of sight, such that the reddened stars on the far side of 
the gas layer have the same intrinsic GMD as the unred¬ 
dened stars; and (2) that the distribution of extinctions 
is log-normal, we use the unreddened GMD and eqn. [I]to 
generate models for the NIR GMD for arbitrary values of 
the median extinction Ay, the dimensionless width a of 
the log-normal distribution of ex tinct ions, and the frac¬ 
tion of reddened stars fred (Sec. 4.3). The first two of 
these parameters characterize the properties of the dust 
distribution, and the third constrains the relative geom¬ 
etry of the gas and stars. We then use a Monte Garlo 
Markov-Chain (MCMC) to evaluate the relative poste¬ 
rior probability that a given set of model para mete rs are 
an acceptable fit to the observed GMD (Sec. 4.4). Fi¬ 
nally, we repeat this process in a series of interleaved, 
oversampled pixels to generate a map of the three pa¬ 
rameters that describe the dust distribution, and their 
associated uncertainties. We now describe the details of 
all the stages in this process, for the specific data in the 
PHAT survey. 

4.1. NIR Data 

We select stars with photometry in the WFC3/IR 
FI low and FI60W filters f rom the simultaneous six- 
filter photometry described in Williams et al. (|20I4). We 
analyze individual “bricks”, which are large contiguous 
rectangular regions consisting of 3x6 WFC3/IR point¬ 
ings that have been astrometrically aligned, merged in 
overlap regions, and culled at the outer edges such that 
each star appears only once in a “brickwide catalog”. Al¬ 
though we include Brick I in fits involving the local stel¬ 
lar surface density to improve continuity, we do not an¬ 
alyze its reddening distribution; this brick is dominated 
by M3I’s high surface brightness bulge and the majority 
of its photometry is quite shallow due to high crowd¬ 
ing. The stellar population in the bulge also contains a 
significant high metallicity population, leading to a very 
broad RGB, limiting our ability to apply the extinction 
mapping technique. Although w e cannot map the dust 
in the bulge using resolved stars, Dong et al. (2014) has 


analyzed its dust distribution using PHAT’s integrated 
photometry in tandem with GALEX and SWIFT obser¬ 
vations. 

We apply further culls to the publicly released pho¬ 
tometry to reduce the number of stars with obviously 
bad photometry (largely due to unresolved blends and 
diffraction spikes from bright stars). While the number 
of such stars is small, their errors typically place them 
in the region redward of the RGB, close to the detection 
limit in FIIOW where they can be mistaken as highly 
reddened stars. Thus, even a small amount of contami¬ 
nation can bias the reddening distribution to spuriously 
high values of the extinction. 

To reduce rare instances of corrupted or highly u ncer¬ 
tain photo metry, we adopt cuts on DOLPHOT’s (Dol¬ 
phin pO^) photometric quality parameters to ensure: 


(I) high signal-to-noise; (2) low probability of signifi¬ 
cant errors in subtraction of neighboring stars; and (3) 
efficient rejection of cosmic rays and of false-detections 
on the diffraction spikes of bright stars. Specifically, 
we require that the signal-to-noise ratio of each star 
was greater than 5 in both NIR filters, that the crowd¬ 
ing parameter (which indicates the amount of contam¬ 
ination from neighboring starp^ and “roundness” pa¬ 


rameter were small (crowdi^ -f- crowd 2 ^ < 2.5^ and 
roundi^ + round 2 ^ < 4.0^). We also restrict the sharp¬ 
ness parameter, which evaluates whether the stellar pro¬ 
file is sharper or flatter than the point spread function. 
We fit an ellipse to the 2-dimensional distribution of 
sharpness in both filters and stars whose distribution in 
sharpi versus sharp 2 fell within an ellipse with major 
axis length 0.35, position angle —128 deg, and axis ra¬ 
tio 0.55. These combinations eliminate clear outliers in 
the space of photometric quality parameters while elimi¬ 
nating very few stars with clearly good photometry (i.e., 
those that fell on the GMD features occupied by the vast 
majority of stars). 

We have considered whether these cuts are unfairly re¬ 
moving reddened stars. We find that there is no a priori 
reason why the reddened stars should be any more af¬ 
fected by the cuts in crowding and sharpness than the 
unreddened foreground, given that the cuts remove stars 
with known photometric limitations that are due to fea¬ 
tures that are unaffected by reddening (i.e., the random 
probability of blending with a star of a small projected 
separation, or of lying on a cosmic ray or diffraction 
spike). We also see no anti-correlation between the num¬ 
ber density of stars and the local extinction. Instead, 
the surface density of stars is very smooth throughout 
the survey area (left panel of Figure where stars is 
defined to be the number density of stars in arcsec“^ 
with 18.5 > FI60W > 21.0 and FIIOW - FI60W > 0.3, 
generated in 15" wide pixels). 

To minimize the effect of incompleteness and magni¬ 
tude biases, we impose additional magnitude cuts beyond 
those implicit in the signal-to-noise cut. We use artifi¬ 
cial star tests to define the 50% completeness limit, mso, 
in each filter (i.e., the magnitude for which a star has a 
50% chance of being detected) for a large number of fields 
across our survey area. Fainter than the 50% complete¬ 
r‘s See|Dolphin| (|2000| for the specific definition of the photomet¬ 
ric quality parameters used here. 




















Dalcanton et al. 



Fig. 3.— Mapping the stellar surface density T^stars (Sec. |4.1| ). (Left) Logio of fho surface density of NIR-detected RGB stars (in 
arcsec“^) with 18.5 > F160VF > 21.0 and F160VF > 0.3, calculated within 15" wide pixels. Outside of the dense inner bulge, the surface 
density falls smoothly with increasing distance from the center, with no obvious “holes” caused by regions of high dust extinction. For 
reference, we plot solid ellipses to represent the isophotes expected for an inclined disk with a position angle of 38.5°, inclination of 74°, 
and major axis lengths of 0.25°, 0.55°, and 1.05°; compared to these locii, it is clear that the internal stellar density distribution of M31 is 
complex, showing multiple components. (Right) A smooth model of the NIR surface density, generated from the sum of 12 Gaussian radial 
surface density profiles, each with a different position angle and inclination. The model reproduces the observed distribution to within 
~15% at all radii. 


ness magnitude, the detection fraction falls dramatically 
while the photometric bias rises, due to an increase in 
the fractional flux contributed by blends with faint un¬ 
detected stars. For each filter we fit a polynomial to 
77150 as a function of the local stellar density, as deter¬ 
mined by a multicomponent model fit to the log of the 
stellar surface density in bright RGB stars (i.e., right 
panel of Figure |^. The data and polynomial fits are 
shown in Figure^ We then use the local surface den¬ 
sity to define a spatially-variable magnitude cut at 77750 
that we apply to the photometry. The effect of the 50% 
completeness cut can be seen in Figure which shows 
77750 superimposed on CMDs constructed from low ex¬ 
tinction regions at a range of representative local densi¬ 
ties. We make the F160W cut even more stringent (by 
A^50 ,fi 60VF = 0.5 mag) than shown in Figure]^ to re¬ 
move stars near the bottom of the CMD, where even 
small amounts of reddening would move the star beyond 
the magnitude limit in FllOW; such stars are numerous 
but do not add significant leverage in measuring the red¬ 
dening distribution and thus we wish to minimize their 
impact on the fit to the CMD. 

In addition to culling poorly measured stars, we also 
make an initial correction for spatially-dependent varia¬ 
tions in the photometry across th e WFC3/IR chip. The 
version of the photometry used in Williams et al. (2014) 
does not include spatial variation in the FbF' model in the 
NIR. This can lead to systematic variations in the pho¬ 
tometry across the face of the camera, such that stars 
in one position on the chip will be measured systemat¬ 
ically brighter or fainter than in another. Errors in the 
large scale flat-held can also lead to similar issues. While 
these effects are partially corrected for by DOLPHOT’s 
perturbations of the PSF model, small differences on the 
scale of a hundredth of a magnitude persist. Our sensi¬ 
tivity to dust depends on small shifts in the NIR color, 
and thus failure to correct for these effects will limit the 
sensitivity we have to low levels of attenuation. 

We develop a position-dependent correction map for 



Fig. 4. — The depth of the photometry as a function of radius. 
The blue and red points show the 50% completeness limit mso as a 
function of logio of the local stellar surface density in bright RGB 
stars, for the FllOW and F160W filters, respectively. The solid 
lines show 10th order polynomial fits to the data, used to interpo¬ 
late to an appropriate limiting magnitude at each position in the 
galaxy. The PHAT photometry is crowding-limited, and thus the 
limiting magnitude is a strong function of the local stellar surface 
density. The apparent roll-off in the inner disk (at high surface 
densities in the left side of the plot) is an artifact of the extraordi¬ 
narily high crowding levels in the inner bulge, which lead to biases 
in the photometry; while we include these regions for completeness 
here, we do not actually analyze their Ay distributions. 


the color by analyzing fi elds in the outer ga laxy that 

(2014) emission- 


Draine et al. 


have low extinction on the 
based dust map. We use 54 individual flelds in Bricks 
2, 6, 8, 10, 11, 12, 14, 18, 19, 20, 22, and 23, and cal¬ 
culate each star’s color relative to the RGB locus. We 
bin the stars as a function of position on the NIR chip 
in a 40x40 grid, keeping only bright stars that fall in 
positions where the extinction Ay ^emission predicted by 


Draine et al. (2014) is less than 0.15mag, and then cal- 
culate the median of the color offset. The resulting map 
shows coherent structure, with redder colors found for 
the largest values of T, peaking at the middle of the 
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range in X. Bluer colors are found for the lowest val¬ 
ues of F, peaking for large values of X. The amplitude 
of the variation is small, with peak-to-peak variations of 
roughly TO.01 mag, co mparable to the amplitude seen by 


Williams et ah (2014 ; see their Fig. 24). The structure 
appears robust with respect to different choices of mag¬ 
nitude range, binning scheme, mean vs. median, field 
choice, etc. We find similar patterns when using colors 
derived from measuring the magnitude of the red clump, 
although the resulting maps are less accurate due to the 
red clumps larger intrinsic width. 

We apply this correction to all of the photometry be¬ 
fore beginning analysis, choosing to modify the F160W 
magnitude to produce the required color shift. Although 
the pattern in the map of color shifts undoubtedly reflects 
differences in photometry between FllOW and F160W, 
we cannot deduce the magnitude shifts from the RGB 
color alone. In practice, the effect of shifting magnitudes 
by ±0.01 mag is minimal, since these shifts are signifi¬ 
cantly smaller than the magnitude error of the size of the 
magnitude bins we use when fitting the observed RGB. 


4.2. Generating a Model for the Unreddened RGB 

Building a model of the unreddened RGB that includes 
spatial variations is essential to mapping the extinction. 
Our extinction mapping technique relies on detecting 
stars that are “redder than expected”, and thus spatial 
variations in the expected stellar colors can be misin¬ 
terpreted as reddening if the underlying RGB model is 
incorrect. We therefore must build a model for the unred¬ 
dened RGB and red clump that captures as much of these 
variations as possible. 

Unfortunately, stellar population gradients, density- 
dependent photometric quality, and projection effects 
make the GMD structure of the RGB and red clump 
complex and difficult to model. Generating a theoreti¬ 
cal model for the unreddened RGB would require solving 
for the detailed star formation history across the disk to 
account for gradients in age and metallicity, while also 
building a model for projection effects and photometric 
errors. Such a model would be subject to errors in the 
disk model, the characterization of photometric errors 
throughout the disk, and the star formation history itself, 
as well as any errors in the underlying stellar isochrones 
and atmospheric models, which can be particularly sig¬ 
nificant for cool stars in the NIR. 

Instead of trying to simulate these effects, we adopt an 
empirical approach. We isolate low-extinct ion regions 
throughout the galaxy, and then bin the stars in these 
regions to generate a model GMD. We create these mod¬ 
els as a function of the local stellar surface density, which 
captures variations in both the underlying stellar popu¬ 
lations and in the photometric errors. 

Here we discuss the process used to generate the 
spatially-dependent model of the unreddened RGB, the 
properties of the RGB itself, and the associated uncer¬ 
tainties. 


4.2.1. Spatial Variations in the NIR RGB 

Almost all of the stars in the PHAT NIR photometry 
he on the RGB or in the “red clump” found at the red 
end of the horizontal branch. These stars span a range 
of ages (> 1 Gyr) and metallicities, depending on the ex¬ 
act star formation and chemical enrichment history of 


M31. The colors of these stars depend on their age and 
metallicity, with older and more metal rich stars favor¬ 
ing redder colorj^ The RGB’s mean color and width 
can therefore change as a function of position if the past 
star formation history and enrichment varies spatially. 
Metallicity and age gradients are common in disk galax¬ 
ies, and distinct structures such as M31’s spheroid and 
bar may likewise have star formation histories that differ 
from those in the disk. The net result is that the struc¬ 
ture of the RGB is likely to vary with position within the 
galaxy. While the RGB is sufficiently old that dynamical 
evolution will have smoothed out any sharp features, fail¬ 
ure to account for the anticipated smooth variations will 
lead to errors when trying to model the color distribution 
of RGB stars. 

The spatial variations in RGB and red clump structure 
can also be subtly altered by viewing angle. The RGB 
and red clump are dominated by older stars, many of 
which were born in or dynamically heated into a thicker 
disk. When this disk is inclined, a range of radii will be 
present along a given line of sight. The stars at a given 
projected radius are therefore a blend of the stars from a 
range of radii and scale heights. The degree of this “pro¬ 
jection mixing” will depend on distance from the major 
axis. Stars on the major axis come from a range of scale 
heights and azimuthal angles, but will all be at the same 
radius. Away from the major axis, however, projection 
effects will mix together stars from a range in radii, not 
just scale height. This projection mixing will leave a sig¬ 
nature on the structure of the RGB and red clump as a 
function of distance from the major axis, if the age and 
metallicity vary with radius and if the disk is thick com¬ 
pared to the scale length over which the intrinsic RGB 
properties vary. We currently do not try to capture the 
impact of projection mixing, which requires knowledge 
not just of the radial gradients in the stellar populations, 
but on the amplitude of such gradients with scale height; 
we expect that projection mixing will be second-order 
compared to the radial gradient and photometric errors. 

Finally, in addition to spatial variations due to projec¬ 
tion effects and in the intrinsic properties of the RGB and 
red clump, we expect there to be spatial variations in the 
quality of the photometry. The NIR observations of M31 
are highly crowded, which affects the behavior of photo¬ 
metric errors, biases, and depth. In general, all of these 
quantities are adversely affected when the stellar density 
increases, such that photometric errors and biases are 
larger at a given magnitude , and the 50% comple teness 
limit is much brighter (see Williams et al. 2014). We 
therefore expect the RGB to be broader and shallower in 
regions of the highest stellar density. 

To characterize the unreddend RGB as a function of 
position, we use the stellar surface density T^stars to track 
position, rather than projected radii. We choose to work 
in the space of stellar surface density (rather than ra¬ 
dius), because the local surface density is the dominant 
factor setting the quality and depth of photometry for 
crowding-limited images. This variation in photometric 


M31’s older, more metal poor RGB stars have similar NIR 
colors to its younger, more metal rich population (i.e., the bluer 
colors of younger RGB stars are partially cancelled by their be¬ 
ing more metal rich), and thus the NIR RGB is typically narrow. 
This narrowness allows even modest amounts of reddening to be 
detectable. 
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Major Axis Length (degrees) 


Fig. 5. — The mapping between Logio of the surface density of 
bright RGB st ars a nd the distance along the major axis (see dis¬ 
cussion in Sec. |4.2|> . Points are color coded by the standard devia¬ 
tion of the upper RGB, relative to a fiducial locus, as in Figure 
The color coding indicates the position of the major star-forming 
ring and the more minor outer star-forming arm, at ~ 0.75° and 
1.15° from the center, respectively. The outer disk declines 
approximately as an exponential disk, but the inner disk is mul¬ 
tivalued, because the single inclined disk model used to calculate 
the major axis length is not appropriate for the complex inner disk. 
The falloff towards the center (<0.2°) is due solely to increasing 
incompleteness in the highly crowded inner bulge. 


quality has far more of an effect on the CMD morphology 
of the unreddened RGB than do the radial gradients in 
the underlying stellar population (which are small, as we 
show below). Moreover, the local density is an equally 
good proxy for “distance from the center of M31”, allow¬ 
ing it to track stellar population gradients without need¬ 
ing to be tied to a specific model for the 3-dimensional 
stellar distribution; a single disk model is not appropri¬ 
ate where the structure of the galaxy is far more complex 
(e.g., the inner bulg e and bar (e.g., [Bea ton et al.||2QQ7) 


and the outer warp (Innanen et al.||l982|)). 'The decrease 
in the stellar density with radius, and the structurally 
complex inner region can be seen in Figure where we 
plot the log of the surface density in bright RGB stars 
(from Figureas a function of the projected major axis 
for a fiducial inclined disk mode0 


4.2.2. Identifying Candidate Low Extinetion Regions 

We isolate low-extinct ion regions by identifying areas 
where the RGB is narrow and unreddened. We first 
divide the stars within each PHAT “brick” into bins 
15" on a side (^ 57 pc) and then make a broad cut in 
color and magnitude to isolate stars on the upper RGB 
(19.0 < F160W < 22). These regions are larger than our 
fiducial 25 pc analysis regions to ensure there are enough 
stars on the upper RGB to provide reliable measurements 
of the RGB’s width and color. We then calculate mean 
color of the RGB relative to a fiducial RGB locus (ap¬ 
proximated as FI low — F160W = uq + ai(F160W — 

26 yyg refer to the projected major axis length as “radius” 
throughout this paper. We adopt a fiducial inclined disk when con¬ 
verting position into radius, assuming PA=38.5°, inclination=74°, 
and center (q;o,(5o) = (10.6847929°, 41.2690650°) in J2000, which 
approximates the disk’s isodensity contours in the Spitzer 3.6p,m 
image of M31. However, M31 has significant structural complex¬ 
ity, and thus this projected radius is not necessarily an accurate 
measure of the distance to the center of M31. 


ino) + a2(F160W — mo)^ where mo = 22.0, ao = 0.752, 
ai = —0.1, a 2 = —0.006). We also calculate the standard 
deviation of the distribution of color differences between 
RGB stars and the fiducial locus; we refer to this stan¬ 
dard deviation as the “width” of the RGB in subsequent 
plots. This definition of “width” is somewhat sensitive to 
differences in slope between the fiducial and the actual 
RGB. It therefore does not solely indicate broadening, 
given that an RGB with a different slope than the fidu¬ 
cial could also have a significant width by this measure, 
even if it were intrinsically narrow. However, such slope 
differences are modest and vary smoothly with radius 
(Fig.left panel), and thus the “width” will still iden¬ 
tify aTocally broader RGB. 

Figure shows the resulting maps of the width of the 
RGB and the mean color shift of the stars with respect to 
the fiducial RGB locus. The main gas-rich star-forming 
rings clearly show up as regions that have larger RGB 
widths and redder colors than is typical for their location 
in the galaxy. The shift in mean color is also largest 
eastward of the major axis, which is the near side of 
M31. In this region, a much larger fraction of stars are 
reddened, as is obvious from optical images of the galaxy. 

We isolated the low-extinction regions from these maps 
as follows. First, we assembled a list of candidate low- 
extinction regions for each brick. Within each brick, we 
tagged regions by their projected major axis length for 
the fiducial inclined disk model, and sorted these sub- 
regions into 20 bins of radius within each brick. Then, 
within each of those radial bins, we identify candidate low 
Ay regions by tagging the 20% of subregions that have 
the narrowest RGBs, as characterized by the standard 
deviation in the color of their RGB relative to the fiducial 
RGB locus. This process guarantees that candidate re¬ 
gions are identified uniformly throughout the brick. We 
repeat this process for every brick and then merge the 
candidate low-extinct ion regions into a list that is guar¬ 
anteed to sample the entire area covered by the PHAT 
data. 

4.2.3. Produeing an RGB Model as a Funetion of Loeal 
Surfaee Density 

After merging the candidate low-extinct ion regions 
from all bricks, we analyze the resulting list to find the 
lowest exinction subregions associated with a given stel¬ 
lar surface density. 

From the list of low extinction candidates, we select 
the regions that have the narrowest RGBs as a func¬ 
tion of local surface density, by including any region that 
falls below the heavy black lines in the left panel of Fig- 
ure {ctrgb = 0.071 + 0.0125(logio + 0.2) and 

ctrgb = 0.071 + 0.0235(logio 'Egtars + 0.2), for the high 
and low density regions, respectively). We then impose 
an additional cut on the color of these regions by itera¬ 
tively fitting a fourth order polynomial to the mean RGB 
color relative to the RGB fiducial, and selecting the sub¬ 
set of regions with the narrowest RGB that are also no 
more than 0.5 standard deviations redder or 7 standard 
deviations bluer than the polynomial fiQ the selected 

The width of this selection region is ~ 0.04 magnitudes in 
color, which is larger than the systematic magnitude variation 
that we see across the face of the WFG3/IR chips (right panel 
of Figure]^; correcting the WFG3/IR calibration would allow us 
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Fig. 6. — Maps of the width and color of the upper RGB (Sec. |4.2.2|. (Left) Standard deviation of bright (19.0 < F160VF < 22) stars 
around a fiducial RGB locus, as a function of position. The RGB is broader in regions of higher extinction. The typical width of the RGB 
also increases towards the inner galaxy due to larger photometric errors from stellar crowding and undetected blends, and potentially due 
to a broader range of intrinsic properties in the underlying stellar population. (Right) The mean color of bright (19.0 < F160VF < 22) stars 
relative to a fiducial RGB locus, as a function of position. The mean colors are redder in regions of high extinction, as expected. There 
are also radial trends towards bluer relative colors in the inner disk. However, this does not indicate a systematic bluing of the population, 
and instead refiects increasing mismatch between the RGB and the fiducial locus. Note also that there is some repeating structure in the 
mean RGB color that mimics the positions of the WFG3/IR fields. This structure is likely to be due to residual < 0.015 mag errors in the 
WFG3/IR fiat fields, or in the spatially-varying PSF model. Ellipses are at the same location as in Figures]^ 


color range is indicated with light solid lines in the right 
panel of Figure The red selection limit culls any po¬ 
tentially reddened regions and the blue selection limit 
reduces fields with large contamination on the blue side 
of the RGB from asymptotic giant branch or red core 
Helium burning (RHeB) stars. Finally, we compare the 
surviving regions t o the dust map recently published by 
Draine et al. ( 2Q14[ ), and reject any regions that fall where 


the implied extinction is greater than Ay > 0.175. The 
stars that fall in the remaining regions are then tagged 
with the local stellar surface density, and then used to 
generate model unreddened CMDs. 

The regions that survive the cuts are highlighted in 
dark blue in the right panel of Figure There are 
^ 153,000 stars total in these regions, sampling the 
full range in local surface density, although not uni¬ 
formly; for example, there are no low extinction lines 
of sight identified in the dusty lOkpc star forming ring 
(—0.5 < logio ^ stars ^ ~0-2), and there are many more 
low extinction regions that sample the diffuse outer disk. 
The spatial distribution of these regions is shown in Fig¬ 
ure superimposed on a map of the width of the upper 
RGB. By design, these regions li e in the lowest extin ction 
regions of the recently published Draine et al. (2014) dust 
mass map, derived from models of dust eniission. We 
note, however, that our selection does not ensure that 
these regions truly have zero extinction. The RGB has a 
finite width due to photometric uncertainties and stellar 
population effects, and we cannot detect extinction vari¬ 
ations that broaden the RGB by amounts smaller than 
the intrinsic width. 

To generate the final model GMD for unreddened stars 
as a function of local surface density, we group the unred¬ 
dened stars into bins of stars- We first rank the 


to narrow this selection region further, potentially improving the 
measurement of Ay at low extinctions if the true population width 
is in fact narrower. 


regions by their local surface density, and then use a 
sliding, adjustable bin to generate CMDs which each 
contain at least 2500 bright RGB stars in the range 
19 < F160W < 21.5). Each of the resulting set of bins 
contain comparable numbers of stars, but do not sam¬ 
ple equal ranges of local surface density; in the low sur¬ 
face density outer disk, one must merge stars in a wider 
range of log^^g ^stars lo reach the same total number 
of stars. We also allow sequential bins to overlap by 
one third, to provide a smoother interpolation over local 
density. To avoid having unnecessarily high numbers of 
nearly identical adjacent bins, we then merge neighbor¬ 
ing density bins together until there is a step of at least 
Alogi^T^stars = 0-05 between each bin. We then keep 
a record of the mean and the range of log^^g ^stars asso¬ 
ciated with each bin. When fitting data, we adopt the 
model GMD from the bin with the closest mean surface 
density to that of the region being analyzed. 


4.2.4. Properties of the RGB in Low Extinetion Regions 

Four representative CMDs of the low extinction re¬ 
gions, in narrow ranges of local surface density, are shown 
in Figure These CMDs are extremely clean, showing 
only a narrow RGB sequence. The CMDs also demon¬ 
strate the quality of the photometric cuts, as measured 
by the nearly total absence of spurious detections near 
the magnitude limit. 

In Figure prl we show the properties of the unreddened 
RGB as a function of magnitude (calculated in bins of 
0.15 mag), in each of the resulting bins of local surface 
density. The mean color of the RGB (left panel) is sys¬ 
tematically bluer and the RGB slope is steeper towards 
low local surface densities (i.e., towards the outer disk), 
as would be cha racteristic of an increas ingly metal poor 
population (see Gregersen et al. 2015). However, the 
shift in RGB morphology is extremely small between ad¬ 
jacent bins, such that our choice to bin in log^g ^stars 
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Fig. 7. — Identifying low exinction regions using the width and mean color of the RGB, 
deviation of the upper RGB, relative to a fiducial RGB locus, as a function of the log of the local stellar surface density of bright 
stars, color coded by the mean color of the RGB relative to the same fiducial RGB locus. Solid lines show the selection limit below which 
we identify the lowest extinction regions at each stellar surface density. (Right) The mean color of the RGB relative to the fiducial RGB 
locus, as a function of the log of the local stellar surface density of bright RGB stars, color-coded by the width of the upper RGB (defined 
as the standard deviation, relative to the same fiducial RGB locus). Light solid lines indicate the range in color used to exclude regions with 
possible extinction (upper curve) or high contamination from RHeB and/or AGB stars (lower curve) from the select ion region identified 
in the left panel. Heavy solid points indicate regions that survive both ciRs, and which also have low extinction in the|Draine et al. (|2014|) 
dust map. The location of these points on the disk are shown in Figure 



Fig. 8. — Location of the lowest extinction regions (blue dots), 
superimposed on a map o f the width of the RGB (e.g., left panel 
of Figur_£^ se^Sec. |4.2.2|). Ellipses are at the same location as in 
Figures [3 and 

leaves no detectable imprint on the extinction maps that 
are eventually derived from the unreddened CMDs. 

The right panel of Figure shows the standard de¬ 
viation of the color of the RGB as a function of magni¬ 
tude. The RGB becomes increasingly narrow at bright 
magnitudes and towards the outer disk, as would be ex¬ 
pected for decreasing photometric nncertaintiej^ and a 
narrower mixture of stellar populations. There is also 
a clear floor at ^ 0.025 mag, which is due to a combi¬ 
nation of the instrinsic width of the RGB, photomet¬ 
ric uncertaintes, and systematic calibration errors in the 
WFG3/IR chip (±0.01 mag). The minimum width cor- 


The only exception is right near the RGB tip at F160W ~ 
18.5, where the decreasing effective temperature of the stars 
leads to strong metal-dependent effects on the stellar photosphere, 
broadening the distribution of colors at the coolest temperatures. 


responds to roughly 0.3 mag of extinction, setting a 
threshold below which we are significantly less sensi¬ 
tive to broadening due to reddening. We could poten¬ 
tially improve our sensitivity to low-but-not-zero extinc¬ 
tion with improved WFG3/IR calibrations if systematics 
dominate, but may not be able to if this floor is the in¬ 
trinsic width of the underlying population. 


4.2.5. Limitations in Building an Unreddened RGB 

There are a number of unavoidable limitations in our 
procedure for generating a model of the unreddened dis¬ 
tribution. While these effects are quite small in terms of 
their effect on the structure of the RGB, they can poten¬ 
tially affect the inferred extinctions. In addition to their 
impact on the model of the unreddened RGB, these ef¬ 
fects will also impact measurements of the redding at low 
extinction. We now discuss each of these in turn. 


The Finite Width of the RGB 

The width of the RGB (Figure 10) creates a limit on 
the lowest extinction we can reliably identify at a given 
location. If the color change produced by a given Ay is 
much smaller than the typical width of the RGB, then 
small redward shifts of the GMD will not be detectable. 
The finite width of the NIR RGB thus directly affects our 
ability to accurately identify the lowest extinction lines 
of sight. 

Although we have taken great pains to identify the 
lowest reddening regions throughout the galaxy, there is 
no guarantee that these are truly zero reddening lines of 
sight, due to the fact that low level of diffuse dust may 
not produce measurable shifts in the NIR colors and/or 
width of the RGB. We have attenipted t o mitigate for 
this effect by using the Draine et al. (2014) dust emission 
maps to place more stringent limits on the local extinc¬ 
tion than are possible with the NIR GMD alone. How¬ 


ever, due to background subtraction, the [Draine et al. 
(2014) maps have uncertainties at low extinction as well, 
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Fig. 9. — CMDs of stars in 4 u nred dened regions, spanning a range of local RGB surface density Tigtars (Sec. |4.l| and/or deprojected 
semi-major axis distance (r) (Sec. |4.2] >. The black lines indicate the 50% completeness limits adopted at that stellar surface density. Note 
the near absence of stars redward ot the CMD, indicating the quality of the photometric cuts and of the identification of low extinction 
regions. Red arrows show the expected effect of 5 magnitudes of extinction in V. Extinction will cause fainter stars to drop below the 
completeness limit at lower Ay than for brighter stars. The effect of the magnitude cuts on the range of observable Ay is more severe in 
the inner disk, where the cuts are closer to the tip of the RGB. The RGB in this region is also intrinsically wider, making it more difficult 
to measure small offsets due to reddening. 


such that some stars with Ay > 0 may still remain in 
our “zero reddening” model CMD. The resulting model 
CMD would be biased slightly to the red, and would 
lead to underestimating Ay when modeling the observed 
CMD. 

We can estimate the amount of reddening that could 
be present in our “unreddened” CMD by inspecting the 
widths of the RGB at bright magnitudes, shown in the 
right panel of Figure If the width of the RGB is 
o'Fiiow-FiGOW ^ Q; Q3 over a reasonably well populated 
part of the RG^^ then one could expect to identify 
extinction at the level of Ay = 0.25, indicated by the 
vertical arrow. We therefore expect that the model of 
the unreddened CMD could be “contaminated” with red¬ 
dened regions with no more than this level of extinction. 
Inspection of the RGB widths in Figure [^suggests that 
identifying low reddening regions will be harder in the in¬ 
ner galaxy, where the RGB appears to be wider, both due 

This plot also suggests that the well populated red clump at 
F160W ^ 22.6 is of somewhat limited utility for detecting small 
extinctions, because the intrinsic width of the RGB+red clump 
feature is much broader than the RGB alone. 


to higher photometric errors from crowding, and from the 
wider range of metallicities present as one nears M31’s 
bulge. However, the overall le vel of extinction in the in¬ 
ner galaxy appears to be low (jDraine et al.||2014|), which 
may perhaps mitigate against this ettect. 

Including bluer filters in our analysis could potentially 
improve our ability to correctly identify low extinction 
regions, while also giving more sensitivity in measuring 
Ay. In Appendixwe discuss the trade-off between the 
intrinsic width of me RGB and the sensitivity of RGB 
color to extinction in other PHAT filter combinations. 


Stellar Population Gradients 

Spatial variation in M31’s underlying stellar popula¬ 
tion pose another possible obstacle to creating a realis¬ 
tic model of the unreddened GMD. We construct mod¬ 
els from stars in individual low reddening regions, but 
frequently, these are far from the individual pixel be¬ 
ing analyzed. Thus, if there are spatial variations in 
the stellar population, we may sometimes use a model 
GMD that does not correctly reflect the local popula¬ 
tion. To first order, this should not present a severe 
problem because we match the local stellar density of 
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Fig. 10.— The morphology and width of the unreddened RGB as a function of position within M31, as characterized by the stellar surface 
density T^stars ■ (Left) Color locus of the unreddened RGB as a function of magnitude. Each line corresponds to a bin of log^^Q Tigtars (in 
units of arcsec“^), which can be mapped to major-axis length using Figure 1^. The mean color was calculated in bins of 0.15 mag, above 
the limiting 50% completeness magnitude shown in Figure ^ The morphology of the RGB varies slowly and smoothly from the inner 
regions (pink) to the outer disk (dark blue). The RGB can K traced to fainter magnitudes in the less crowded outer regions, due to the 
smaller photometric errors and fainter crowding limit. (Right) Standard deviation of stars around the mean color locus as a function of 
magnitu de, ca lculated in the same bins as in the left panel. The arrow indicates the expected shift in color for an extinction of Ay = 0.25. 
See Sec. |4.2.4] for a discussion of this figure. 


the model to the pixel being analyzed, which isolates 
regions at comparable deprojected radii, and thus of 
similar underlying stellar population. Empirically, Fig¬ 
ure ^ shows that the RGB color and width is a very 
weakmnction of local stellar density, particularly where 
logio < —0.2 arcscc”^. Thus, systematic offsets in 

the structure of the RGB are likely to involve color dif¬ 
ferences of only a few hundredths of a magnitude, which 
would produce a very small bias in the modeled RGB. 
The only place where we anticipate that stellar popula¬ 
tion gradients might significantly affect the model unred¬ 
dened GMD is in the inner regions of the galaxy where 
the stellar populations change more rapidly with radius, 
and where the structure of the galaxy is complex (due to 
M31’s bar, and multi-c omponent bulge; see Figureand 
Gregersen et al.| (|2015|). 

A potentially more significant impact of stellar popula¬ 
tion gradients comes through their effects when viewing 
an inclined disk (as in M31). Projection effects lead to 
pixels along the minor axis having stars from a wider 
range of radii than along the major axis. The net result 
is that while matching of the stellar densities correctly 
matches the photometric properties of the RGB, it may 
not successfully match the underlying color and shape of 
the RGB, due to a mismatch in the stellar populations 
between the model and the pixel being analyzed. 

Foreground Extinction from the Milky Way 


The final limitation when constructing a model for the 
unreddened GMD is that there may be foreground dust 
extinction from the Milky Way or from a diffuse dust 
halo in M31 itself. This foreground would not affect our 
results if it were uniform, because our measurement is 
fundamentally relative; if both the reddened and “unred¬ 
dened” stars are affected equally, the measurement of the 
reddening relative to the unreddened model is still sound 
even if the “unreddened” model has been affected by dust 
in the Milky Way foreground. If the foreground dust is 


not uniform (which indeed is most likely), then it broad¬ 
ens the RGB in the model for the unreddened GMD, due 
to including lines of sight with different reddenings. The 
net result will be to reduce sensitivity to low levels of red¬ 
dening within M31. Unfort unately, one cannot use maps 
of the Milky Way du st (e.g., Schlegel et al.IlQQS Schlafiy 


& Finkbeiner 2011) to correct tor the foreground dust 
since these maps fail when there is a bright background 
source of FIR emission (such as M31) that dominates the 
emission from the Milky Way. 

While the list above describes the principal limitations 
in constructing the unreddened GMD, there are addi¬ 
tional uncertainties associated with th e mod el fitting it¬ 
self. We discuss these below in Sec. 5.3.2[ and make 
empirical tests of systematics in our final maps in Sec. 


4.3. Generating the Model Reddened CMD 

Our data is set of F160W magnitudes rui and 
FllOW—F160W colors q. We therefore need to trans¬ 
late our adopted model for the expected distribution of 
extinctions pa and/or reddenings pg into a model of 
the expected color magnitude diagram for a given set 
of model parameters 0 = {Ay, cr, fredi given the em¬ 
pirical model for the distributions of colors and magni- 
tnrlGs of unreddened RGB stars that we derived in Sec. 
4.2.3 In other words, we wish to derive Pcmd(c, m|^), 
the probability of finding a star at color c and magnitude 
m, given 0. This model will depend on the local surface 
density stars ^ but we do not make this surface density 
dependence explicit, for notational compactness. 

Before preceding with the calculation of pcmd^ we 
take two steps to increase the computational efficiency. 
First, we bin both the data and the models into bins of 
color and magnitude. We use a grid spacing of 0.015 mag 
in color and 0.2 mag in magnitude. This choice guaran¬ 
tees that there are at least 2 pixels across the la width 
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of the narrowest RGB, effectively Nyquist sampling the 
width. Experiments with smaller bins produced no no¬ 
ticeable changes in the resulting extinction maps, while 
dramatically increasing the computation time. 

Second, we translate all F160W magnitudes into 
“reddening-free” magnitudes defined as 

q = F160W - ((FllOW - F160W) - cq) 

^_ Apieow / Ay _ (9) 

{Afiiow / Ay) — {Apimw / Ay) 

where Cq is an arbitrary color for which q = F160W. 
With this transformation, increases in extinction will 
cause a star’s color to redden, while leaving the 
reddening-free magnitude unchanged. With this skewed 
version of the CMD, changes in the extinction move stars 
horizontally in the grid, rather than diagonally (i.e., the 
reddening vector becomes horizontal). We can then do 1- 
dimensional convolutions when calculating the reddened 
CMD, rather than more computationally expensive 2- 
dimensiona]^ convolutions. We will therefore calculate 
PCMD{c,q\0) rather than pc'md(c, m|<9). 

The next step in generating pcMD is to define what 
region of the CMD will be fit. We restrict our analy¬ 
sis to regions occupied by RGB stars by generating a 
mask. The blue boundary of the mask is defined to be 
2 .5cr blueward of the RGB (Fig. [Iq|) to reduce contribu¬ 
tions from AGB and RHeB stars. Phe upper boundary 
is fixed at g = 18.5, which is approximately the tip of the 
RGB. The lower faint boundary and the right diagonal 
boundaries were set by the adopted F160W magnitude 
limit (m5o,Fi60W~0-5) and the 50% completeness limit in 
FllOW, respectively, both translated into the reddening- 
free CMD. All normalizations and comparisons with data 
are restricted to be within the unmasked regions. 

Within each bin of local stellar surface density, we grid 
the unreddened stars into a CMD where the reddening- 
free magnitude q has replaced F160W, eliminating stars 
that fall in the masked region. We then normalize the 
binned CMD so that it integrate to one over the un¬ 
masked area. We refer to the final binned, masked, nor¬ 
malized probability distribution for unreddened stars as 

Punreddenedi,^’) Q)‘ 

We then generated the probability distribution for 

the reddened stars, Preddened{c^ Q\Ay ^ a)^ by convolving 
the binned unreddened CMD with a 1-dimensional log¬ 
normal kernel, where the properties of the kernel are 
set by Ay and a. We then added the reddened and 
unreddened model CMDs together to generate the model 

for the combined CMD Punreddened+reddened{c, q\0) after 

weighting them by appropriate fraction of reddened stars: 


Punreddened-\-reddened{_^-) Q\^) — (1 fred)Punreddened{^’)Q) 

T fred Preddened (^5 Q. \ Ay , (j) • 

( 10 ) 

We then reapply the mask to Punreddened-\-Teddened S'lld 
renormalize. 

We also include a third component to model the noise 
from potentially bad photometry. Occasionally, spurious 
sources will be detected with colors and magnitudes con¬ 
sistent with being reddened RGB stars. These sources 


are rare, but not impossible, and thus we need to include 
an additional component in our model that allows there 
to be a small chance of an individual source being spuri¬ 
ous, rather than requiring every star redward of the RGB 
to be due to dust reddening. To model the CMD contri¬ 
bution from noise, Pnoise^ we identify stars that are more 
than 3 . 5(7 redward of the mean RGB color in regions that 
were identified as being “unreddened”. These red stars 
are unlikely to be due to intervening dust, and instead are 
due to occasional photometric errors when doing crowded 
field photometry. These anomalous red stars are then 
gridded into a “noise CMD”, Pnoise(c, g), and then re¬ 
added into the sum of the unreddened+reddened models 
in the proportion with which they were found in the orig¬ 
inal unreddened CMD. Letting fnoise be the fraction of 
stars in the catalog of unreddened stars were classified as 
noise. 


Pcmd{c, q\0) = (1 - /n 


s) Punreddened+reddenedip'i Q\^) 
T f noise Pnoise (c, q) 

( 11 ) 


The fraction of pixels in the noise model was never more 
than 1.5% in any bin of local stellar surface density, and 
was less than 0.5% in 3/4 of the bins. 

Finally, we introduce a fourth parameter Sc that allows 
the color of the unreddened CMD to shift by a few hun¬ 
dredths of a magnitude, with the goal of absorbing any 
residual issues with the spatially-dependent WFC 3/IR 
calibration/photometry correction discussed in Sec. 4.1 
Bec ause the photometric correction we applied in Sec. 
|4.1| was calculated on a coarse grid of chip position, it 
does not accurately trace rapid changes with position, 
particularly those near the edges of the chip; these resid¬ 
ual issues can be seen as slight bands in Figure that 
trace the edges of the NIR chip in adjacent observations. 
In addition, a handful of fields show small global shifts in 
RGB color, most likely due to a small error in the aper¬ 
ture correction that was applied to the entire field. Both 
of these residual systematic photometry errors translate 
directly into features on the resulting extinction map if 
not corrected for. 

An example of an unreddened and reddened CMD are 
shown in Figure El 


4.4. Fitting the Reddened CMD 


The model CMD shown in Figure 11 has three ad¬ 


justable parameters — the median Ay of the log-normal 
distribution of extinctions, the dimensionless width a of 
the same distribution, and the fraction of reddened stars 
fred (eqns. 0 & - in addition to one parameter Sc to 

help absorb any residual issues in spatially variable pho¬ 
tometry. We solve for these parameters and their uncer¬ 
tainties within small spatial regions of the PHAT imag¬ 
ing data, using MCMC fitting to find the most probable 
fit and the associated uncertainties, subject to sensible 
prior probability distributions (“priors”) on the values of 
the p arameters. The choice of priors is discussed in Sec. 
14.4. II below. 

Specifically, we use Pcmd{cq\^) derived in the previ¬ 
ous section to calculate the probability of finding a star 
at color c and reddening-free magnitude g, for a given set 
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Fig. 11.— Comparing a binned unreddened CMD (left) to a reddened model including noise (right). Star symbols show the location of 
stars found in a 25 pc pixel in Brick 15. The best fit parameters for that pixel are used in the reddened model. Note that the vertical 
axis is not F160W magnitude, but instead is the “reddening-free” F160W magnitude q (eqn. [^, such that changes in the amount of 
extinction/reddening moves stars horizontally. 


of parameters 0 = {Ay, cr^ fred^ assuming a model for 
an unreddened CMD at the local surface density T^stars^ 
The total probability distribution function ptotai of 
measuring a set of N independent and identically drawn 
observations {ci^qi} is then 


N 

Pdata{{Ci,qi}\0) = YlPCMD{Ci,qi\0). (12) 

i=0 

We do not concern ourselves with the variation in the 
number of stars per CMD point, given that any variations 
due to the loss of highly extincted stars is expected to be 
extremely small, adding little useful information to the 
fit. 

We are interested not in the probability distribution 
function (PDF) of a given set of observations given the 
parameters of the adopted reddening model, but instead 
are interested in the posterior probability distribution 
function for the reddening parameters 6>, given the set of 
observations {q, g^}. We therefore use Bayes theorem to 
derive 


Pparam{^\\_Ci’! Qi}) ^ Pdatai^^i-! Pprior{^) (f 

N 

'xYlPCMD{Ci,qi\0) Pprior{0), (14) 
i=0 


Way (e.g., Kainulainen et al.||2009), and the dimension¬ 
less log-normal width of the prior on cr is 0.3. This choice 
also restricts the value of a to be greater than zero. We 
adopt a narrow, flat-topped prior for since we ex¬ 
pect the photometric systematics to be small (i.e., as the 
same o rder as the large scale systematics corrected for 
in Sec. 4.1), but do not have a strong prior on the exact 
value. We adopt a functional form similar to a gaussian 
with as = 0.035, but with the exponent raised to the 
fourth power, rather than squared; this change produces 
a peaked distribution with a flatter top and faster fall-off 
than a gaussian. 

When setting the prior for /red, we take a more so¬ 
phisticated approach than for Ay and a. The value of 
/red is naturally bounded by zero and one, but rather 
than having a hard limit at each of these values, we in¬ 
stead adopt a prior that increasingly penalizes values of 
/red that approach these limits. We do so by “regular¬ 
izing” /red by switching to a closely related parameter 
that more easily penalizes the reddening fraction becom¬ 
ing zero or 1, while maximizing the prior probability at 
an arbitrary mean value of (/red)- Specifically, we adopt 
a new parameter x such that 


/red 


1 -he^ 


(In (/red))/(ln0.5) 


(15) 


where Ppriori^) encodes any prior information on the val¬ 
ues of the parameters. 

4.4.1. Choice of Prior Probability Distributions 
We chose specific forms for Pprior to limit the parame¬ 
ters 0 to physically plausible values. We chose the prior 
in Ay to require that the extinction be greater to or equal 
to zero. For a, we adopt a wide log-normal distribution 
where the most probable value of cr is 0.3, which is char¬ 
acteristic of extinction distributions seen in the Milky 


With this change of variables, fred = {/red) when x = 0, 
and fred goes to zero or 1 only when x Too. We then 
adopt a prior that has a mean {x) = 0 and that declines 
to both positive and negative x, which maximizes the 
prior probability on fred to have the proper mean, while 
increasingly (but smoothly) penalizing values of fred that 
approach the extremes of its permitted values. 

Eqn. p!5] requires chosing a value for the expected mean 
value oTfred- Naively one might expect (fred) = 0-5, but 
this is generically only true along the major axis of an 
inclined disk. For a thickened stellar disk, projection ef- 
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Fig. 12.— The log of the priors in f^ed ^ function of Ay^ for an assumed mean reddening fraction of (fred) = 0-2 (left) and (fred) = 0-8 
(center). The right hand plot shows the log-normal prior in cr, for which the width of the reddening distribution is assumed to be ind epe nden t 
of extinction. These priors can be compared to the distribution of recovered values shown in the upper two panels of Fignres and [T7| 


fects lead the stars on the near side of the disk to come 
from slightly different radii than the stars on the far side 
of the disk. The number of stars declines approximately 
exponentially with radius, and thus the fraction of red¬ 
dened stars at some projected radius can be much higher 
or lower than 0.5, depending on whether the stars in front 
of the dust layer are coming from inside or outside the 
nominal projected radius. As a result, for a moderately 
thick inclined stellar disk like M31, the fraction of red¬ 
dened stars should vary significantly from the near side 
to the far sidq^ by an amount that depends on the incli¬ 
nation and on the thickness of the stellar disk relative to 
the disk’s exponential scale length. As such, the appro¬ 
priate value of {fred) should vary with position to avoid 
biasing the derived parameters. 

To create a prior that looks as much as possible like 
the actual data, we developed a simple geometric model 
for the reddening fraction of a thick, inclined disk. We 
adjust the properties of this model iteratively, using the 
observed spatial distribution of fred regions where the 
measurement of fred has small reported errors, and thus 
is essentially unaffected by the prior. Between iterations, 
the only significant change in the maps of fred are in re¬ 
gions of very low extinction where fred is largely uncon¬ 
strained by the data, and instead gravitates to the value 
of {fred) set by the prior. The fit converged after only 2 
iterations (verifying our assumption that the measured 
values of fred in low uncertainty regions were unaffected 
by the choice of prior). The resulting model for {fred) 
has a position angle of 37°, inclination of 78°, and a ra¬ 
tio of vertical to horizontal exponential scale heights of 
hz/hr = 0.15. The residuals from this simple model are 
approximately A fred = ±0.1. 

After setting the most probable value of fred as a 
function of position and extinction, we set the width of 
the prior probability distribution by using an asymmet¬ 
ric split normal distribution in x. We keep the priors 
broad, to avoid biases in regions where the simple geo¬ 
metric model for fred is not ideal, and to account for our 
more uncertain knowledge of the filling factor of molecu¬ 
lar clouds f fill = *^cloud/'^pixel') whcrC j^eloud and *^pixel 
are the areas of the gas cloud and the analysis pixel, 
respectively. We narrow the prior for low extinctions 
{Ay < 0.3), where the RGB is only slightly broadened, 
rather than split into distinct peaks, since these regions 
will not have sufficient information to constrain fred reli- 

As indeed has long been evident in optical images of M31. 


ably without the prior. Examples of the resulting priors 
for regions with {fred) =0.2 and 0.8 are shown in Fig. 

m 


4.4.2. Fitting Parameters in a Pixelized Map 


Using these priors and the likelihood (eqn. we use 
a MCMC sampler to characterize the posterior prob- 
ability distribution for 0 (eqn . 13). We use emcee 
(Foreman-Mackey et al. 2013), which implements an 
attine-invariant ensemble sampler to efficiently sample 
the distrib ution using a set of coup led MCMC chains 
(“walkers”; Goodman fc Weare||201Q ). 

We run the MCMC sampler within pixels defined by a 
dense spatial grid covering the survey area. Our fiducial 
grid uses square 6.645" pixels, which correspond to 25 pc 
at the distance of M31 (776 kpc). M31 is highly inclined, 
and thus projection effects may make the effective spa¬ 
tial resolution of the survey significantly worse than 25 pc 
in the direction paralleling the minor axis. However, al¬ 
though this degredation in resolution is generically true 
for disk structures, there is no strong evidence in the 
Milky Way that molecular clouds are flattened in the 
same sense as the global galactic disk, rather than being 
3-dimensional objects with no preferred direction to their 
structure. In this latter case, projection does not actu¬ 
ally affect the physical resolution; as an extreme exam¬ 
ple, if molecular clouds were spherical, they would appear 
to have the exact same size when viewed at any angle. 
We therefore adopt 25 pc as the effective resolution of 
the dust map, in spite of the inclination of M31’s disk, 
but recognize that the exact physical resolution could be 
worse in the direction perpendicular to the major axis. 

The adopted 25 pc grid size balances the accuracy in 0 
(which decreases when the re are few stars per pixel, as 
discussed below in Sec. 5.3) and the spatial resolution of 
the resulting map. We have found that we can produce 
reliable maps at significantly higher resolution, but only 
in regions where a large fraction of the stars are reddened. 

We oversample the 25 pc grid by rerunning the fitting 
procedure at 4 dithered locations, shifted by 0.5 pixel. 
We then interleave the resulting 4 dithered maps, pro¬ 
ducing pixels of 12.5 pc on a side. This sub-sampling is 
roughly analogous to Nyquist sampling a map with 25 pc 
resolution. Each of the interleaved pixels shares 50% of 
the stars with the adjacent pixels with which it shares a 
border, and 25% of the stars with the pixels with which it 
shares a corner. Thus, pixels are not independent of their 
immediate neighbors, but are completely independent of 
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all other pixels. 

Within each 25 pc spatial pixel, we use 50 sets of cou¬ 
pled MCMC chains, and take 150 steps to localize the 
chains near the parameter values that maximize the pos¬ 
terior probability. We then run the sampler for a final 
15 steps (i.e., 50 x 15 samples) to sample the posterior 
distribution. To deal with memory limitations and allow 
parallelization, we analyze the pixels in batches, defined 
by the area of individual survey bricks; when a pixel over¬ 
laps adjacent bricks, we use the results for the pixel that 
contained the largest number of stars. Running the sam- 
per for all 22 bricks at all 4 dither positions takes over a 
week when running on 50 cores. 

A representative example of the resulting distribution 
is shown in Figure for the same model shown in Fig¬ 
ure [T^ The distribution is smooth and unimodal, as 
has been the case for all distributions that we have in¬ 
spected individually. Because of the simplicity of these 
distributions, we do not save the full probability distribu¬ 
tion at each grid point, but instead record the values of 
the parameters that maximize the posterior probability. 
We also marginalize the distributions for each parame¬ 
ter, and record the values corresponding to 16%, 50%, 
and 84% of the parameter distribution (i.e., the median 
and the Tier points for a Gaussian distribution). When 
quoting a single value for the uncertainty in each param¬ 
eter (rather than a range), we define the uncertainty in 
the Tth parameter to be /S.0i = (6>i^84 — ^ 2 ,i 6 )/ 2 . 

We also record the equivalent percentile points for 
quantities that we derive from the best-fit parameters 
during subsequent analysis. For example, the mean Ay 
depends on both Ay and a (eqn. [^ , so we calculate the 
marginalized distribution of {Ay) = Ayexp(cr^/2) and 
record the resulting 16%, 50%, and 84% percentile val¬ 
ues, rather than relying on formulas that assume prop¬ 
agation of Gaussian errors in Ay and a. The same de¬ 
terministic translation of uncertainties is used when con¬ 
verting from the regularized variable x back to the red¬ 
dening fraction fred (eqn. 15). 


5. PROPERTIES OF THE REDDENING PARAMETERS 

Before presenting final maps, we use the results of fit¬ 
ting individual PHAT bricks to demonstrate the overall 
accuracy of the method. We show results for two inde¬ 
pendent regions, and discu ss the model param eters ’ spa¬ 
tial distribution (S ec. |5.1| ), correlations (Sec. 5.2), and 
uncertainties (Sec. |5.3[r We then discuss the method’s 


overall susceptibility to systematic errors (Sec. 5.3.2), a 
point we re turn to when comparing to other dust tracers 
in Sec. 16.41 below. 


5.1. The Spatial Distribution of Derived Reddening 
Parameters 


The left panels of Figures [14] and 15 show the spatial 

distribution of Ay, a, and fred for two of the PHAT 
bricks that sample regions with different star formation 
intensities (Bricks 16 and 15, the latter of which con¬ 
tains some of the most intense star formation in the 
PHAT survey area) . Each brick covers an a rea of roughly 
1.5 kpc X 3kpc; see Dalcanton et al. (2012) for details. 

Although the reddening parameters are derived inde¬ 
pendently for each spatial pixel, there are clear coherent 
features in all the derived parameters across the area. 


on scales larger than the one pixel coherence expected 
solely from oversampling the 25 pc pixel grid. The distri¬ 
bution of dust extinction shows clear filamentary struc¬ 
ture throughout the kiloparsec-scale regions, and is rem- 
iniscent of large scale maps within the Milky Way (e.g 


Schlegel et al.|1998l|Froebrich et al.|2007l|Kohyama et al 


2Q13[ [Lombardi et'al. | |2Q11[ ISchlahy fc Finkbeiner | |2Qll 

Nidever et al.|2012p. There is also obvious large scale co¬ 

herence in the spatial distribution of the reddening frac¬ 
tion fred- The width of the reddening distribution a also 
is spatially coherent and follows the structure in the red¬ 
dening distribution, but with a much smaller dynamic 
range. 

5.2. Correlations Among the Reddening Parameters 

Fignres pT| and pT| plot correlations among the different 
reddening parameters, for the maps in Figuresand 
respectively. 

The plots of fred a function of Ay (upper left) reveal 
a number of features that are due to either the dust dis¬ 
tribution or to the fitting process, or both. First, at high 
extinction {Ay > Imag), the mean value of fred shows 
no obvious trend with increasing extinction; this behav¬ 
ior is most noticeable in Figure [T7| where the extinction 
and fred are highest. The simplest interpretation of this 
behavior is that dusty gas largely fills the entire 25 pc 
pixel when the extinction is high, such that the value 
of fred reflects only geometrical effects (i.e., the relative 
position of the stars and the gas along the line of sight). 
In this high extinction regime, the scatter in fre d is due 
both to the uncertainty in measuring fred (Sec. 5.3 be¬ 
low) and to the intrinsic variation in the relative position 
of the gas with respect to the stars. 

At lower extinctions (0.3 < Ay < Imag), the typical 
value of fred tends to fall with decreasing extinction (al¬ 
beit with large scatter, particularly in Brick 15, where 
the value of fred varies significantly across the brick). It 
is unlikely that the relative placement of gas and stars 
along the line of sight would “know” about the overall ex¬ 
tinction of the gas. Instead, this fall-off may reflect that 
dusty gas has an increasingly smaller filling factor as the 
median extinction decreases (see, for example. Figure]^. 
In this regime fred equals the product of the “geometric” 
value of fred,geom that ouc would mcasurc if the gas filled 
the entire pixel and of the filling factor ffm. Thus, if the 
area covered by gas is smaller for lower extinctions, then 
one expects exactly the fall-off we observe. Note that in 
this range of Ay we have not built the decreasing fill¬ 
ing factor into the prior (see Figure 12), and thus the 
decrease effectively reflects the properties of the gas it¬ 
self. On the other hand, some of the reduced scatter 
may be due to the increasing role that the prior plays as 
the information content of the data decreases. At lower 
extinctions, reddening does not cleanly separate the red¬ 
dened and unreddened RGB, making measurements of 
fred more challenging, which increases the weight given 
to the prior and pulls the measured value of fred closer 

fc {fred}- 

At the lowest extinctions {Ay < 0.3mag), the rela¬ 
tionship between fred as a function of Ay changes again. 
In this regime, the width of the NIR RGB is compara¬ 
ble to or less than the broadening expected due to red- 
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Fig. 13. — Samples of the posterior probability distribution function for the width of the log-normal distribution a (left) and the reddening 
fraction (right) as a function of the mean Ay (left) and the median Ay (right). Samples were drawn for the same 25 pc pixel used in 
FigurefTll Bestfit values and the 16%-84% range are fred = 0.84^°'^^, (Ay) = 1.80^q‘24, and a = 0.24^Q'Qg. 


dening. As such, the observed CMD is barely changed 
by the dust, and the constraints on the reddening pa¬ 
rameters becomes extremely weak. As such, the most 
probable value of fred is driven almost entirely by its 
prior. In other words, when the posterior probability 
distribution is very broad, the peak is set primarily by 
the prior probability distribution function, which in this 
case drives fred towards (fred) as the extinction decreases 

below Ay < 0.3 mag. 

Unlike the clear trends between Ay and fred^ any 
trends of a with the other reddening parameters are far 
weaker. Below Ay < 0.75 mag, the distribution of a is 
clearly dominated by the prior. Large values of a are only 
found at higher extinctions {Ay > 1 mag, most preva¬ 
lent in Figure [T7|). Visual inspection of Figures M and 
[T5] confirms that larger values of a are associatedfwith 
regions of high reddening, but Figuresand ^suggest 
this is largely the effect of the fits being able to move 
away from the prior when the signal from reddening is 
strongest. However, while the distribution of a broad¬ 
ens in high extinction regions, it is clearly much more 
likely to scatter to higher values than lower ones, sug¬ 
gesting that the distribution of reddenings can become 
much broader, but rarely becomes narrower, as seen by 
the floor at a ~ 0.25. These high extinction regions with 
broader reddening distributions may reflect differences in 
the underlying turbulent structure, a breakdown in the 
log-normal assumption for the reddening distribution, or 
an increased likelihood of multiple gas layers falling in a 
single analysis pixel when the column density is high. 

5.3. Uncertainties in the Reddening Parameters 

In the right panels of Figures and we plot the 
spatial distribution of the uncertainties in Ay, a, and 
fred (he. Ad = (6^84 — die)/2). The uncertainties are 
clearly not spatially uniform, and instead mimic struc¬ 
tures seen in the distribution of Ay and /red- This cor¬ 
respondence reflects the fact that the accuracy of the de¬ 
rived parameters depends on how cleanly the extincted 
stars separate from the unreddened RGB. Thus, the un¬ 


certainties are lowest when the median extinction Ay is 
high (as in Brick 15; Figure 15). 

The accuracy of the measured parameters also depends 
on the number of reddened stars in the pixel. In Brick 

15 (Figure the uncertainties on the extinction are 
low in the upper right corner, where the fraction of red¬ 
dened stars fred is highest, in spite of the fact that the 
extinction is not particularly high in that region. Brick 

16 (Figure shows a similar gradient in accuracy, but 
at essentially constant fredi here, the fall off in accuracy 
is due to the lower number of stars per pixel overall, 
since the stellar density falls significantly from right to 
left (not shown). Uncertainties are also high for edge pix¬ 
els, which may not be fully covered by the data. These 
pixels disappear when we merge the bricks into a con¬ 
tiguous map, because the edges of the bricks overlap sig¬ 
nificantly. In cases of overlaps, we adopt the fits from 
whichever overlapping pixel had measurements with the 
smallest uncertainties. 

Because the accuracy depends most strongly on the 
extinction of the reddened component, in Figures [Ts] 


and 

and 
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we plot the fractional uncertainties of A Ay/Ay 
Aa/cr, and the absolute uncertainty A fred ^ func¬ 
tion of Ay, for all pixels in Bricks 16 and 15 (whose pa¬ 


rameter distributions were shown in Figures [16] and 17) 

As expected, there is a clear trend towards oetter pre¬ 
cision at higher extinctions. For large Ay, the reddened 
distribution of RGB stars becomes completely distinct 
from the unreddened RGB in the GMD, and the dis¬ 
tribution typically becomes broader, leading to stronger 
constraints on the reddening distribution. At high ex¬ 
tinctions, the typical uncertainties are smaller than 20% 
in Ay, and smaller than ±0.1 in fred- The uncertainty in 
a is comparatively high, being rarely smaller than 20%. 
As with Ay, the fractional uncertainty in a drops some¬ 
what at high extinction, because of the better separation 
of the unreddened and reddened RGB sequences. 

There are rare outliers at high Ay, where the fractional 
uncertainties are much larger than typical. These outliers 
usually correspond to cases where the posterior probabil- 
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Fig. 14. — Left: Maps of the median extinction Ay (top), log-normal width a (middle), and reddened fraction fred (bottom) for Brick 16. 
In regions of low extinction, the values of fred cind a are driven to the most likely value in their priors (= 0.4). Right: Maps of the fractional 

uncertainty in extinction A Ay I Ay (top) and the log-normal width Acr/cr (middle), and the absolute uncertainty in the reddening fraction 
A fred (bottom). Uncertainties are smallest in the regions of highest extinction. In general, the fractional errors in Ay are much smaller 
than in a. 

ity distribution on Ay is very broad due to a very small 
number of reddened stars, as may occur when the total 
number of stars is small, when the reddening fraction is 
very low, and/or when there is a region with bad photom¬ 
etry (typically due to sources detected around saturated 
stars). 

At lower extinctions {Ay ^ Imag), the uncertain¬ 
ties in Ay are more typically 20-30% in regions where 
fred > 0-3, and reach 50% by Ay < 0.3 mag. The typ¬ 


ical uncertainties are larger in regions with even lower 
reddening fractions, due to the smaller number of stars 
in the reddened peak. These regions also have more spu¬ 
rious high Ay points with large uncertainties, because 
the presence of very few red stars can appear to mimic 
high extinction when the prior favors small values of fred- 
The uncertainties in fred and a decrease towards low ex¬ 
tinction, because of the narrowing of the prior rathan 
than because of increasingly constraining data. 


A(j/<7 ^A.y / A.'y 
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Fig. 15.— Same as Figure [l^ but for Brick 15, which falls on the lOkpc star-forming ring on the major axis, and thus samples the highest 
star formation of any of the PHAT bricks. Left: Maps of the median extinction Ay (top), log-normal width a (middle), and reddened 

fraction fred (bottom). Right: Maps of the fractional uncertainty in extinction A Ay I Ay (top) and the log-normal width Acr/cr (middle), 
and the absolute uncertainty in the reddening fraction (bottom). Given that uncertainties are smallest in the regions of highest 

extinction. Brick 15 has low uncertainties overall. 


A comparison of Figures p!8| and also shows the ef¬ 
fect of fred on the uncertainties. Tme reddened fraction 
in Brick 16 (Figure is much lower than the typi¬ 
cal reddening fraction m Brick 15 (Figure 19), leading 
to lower precision at the same median extinction. At 


Ay = 2 mag, the fractional accuracy in the extinction is 
^ 12% in Brick 16, whereas it is better than 6% in Brick 
15, which has a reddening fraction that is nearly twice as 
high (see Figures and|l7|) . This difference reflects the 


larger number of stars in the reddened peak when fred is 
large. 

To better show the impact of the stellar density on the 
uncertainties, the right panels of Figures and show 
A fred and Aaja as a function of the median extinction, 
color coded by the number of stars in the pixel. There 
is a clear tendency for pixels with fewer stars to have 
broader posterior distributions and thus larger uncertain¬ 
ties. For a fixed pixel scale, this trend indicates that the 


A(j/<7 ^A.y / A.'y 
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Fig. 16.— Joint distributions of derived parameters (Ay vs fred: vs cr, and fred vs cr, clockwise from top left) for each pixel in Brick 

16, color-coded by the parameter labeling the figure; L is the likelihood calculated for the fit. At low Ay (< 0.5), the values of fred cind 
(7 are clearly pulled towards the maximum probability of the assigned priors (= 0.4). The bottom right shows fred vs a again, but now 
color-coded by the log likelihood of the model fit. 


uncertainties will be larger in the outer disk, where the 
stellar surface density is low. This trend also suggests 
that better spatial resolution could be achieved in regions 
with high reddening fractions, where we can reduce the 
number of stars per pixel in exchange for increased res¬ 
olution, while still preserving scientifically useful results. 
We have not chosen to adopt spatially-variable resolu¬ 
tion elements here, to allow easier comparisons to global 
emission maps of M31’s ISM. However, we will be pursu¬ 
ing higher resolution analyses of specific high-extinction 
subregions in future work. 

Finally, we stress that the above uncertainties are ran¬ 
dom, and not systematic. However, in the low extinc¬ 
tion regime where systematic biases could potentially 
become important, the random uncertainties are large. 
Thus, while we expect biases at low extinctions (e.g.. Sec. 


4.2.5), we also know not to over-interpret these extinc¬ 


tions because of their significant random uncertainties. 




















Dust Mapping in M31 


23 



Fig. 17.— Same as Figure |l6| but for the high star-formation region covered by Brick 15. Joint distributions of derived parameters (Ay 

vs fredi VS (7, and fred VS (T, clockwisc from top left) for each pixel in Brick 15, color-coded by the unplotted parameter (denoted by 
the legend). The bottom right shows f^ed vs a again, but now color-coded by the log likelihood of the model fit. 
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Fig. 18.— Uncertainties /AAy jAy^ Acr/cr, and Af^ed (clockwise from top left) as a function of median extinction Ay, for each pixel in 
Brick 16, color-coded by the reddening fraction (upper left), n umb er of stars fit per pixel (right panels) and the log likelihood of the model 
fit (lower left). As expected from visual inspection of Figure uncertainties increase strongly at low extinctions, where the reddened 
RGB is not cleanly separated from the unreddened RGB. In regions of modest or high extinction (> Imag), fractional uncertainties are 
typically lower for the median extinction than for the width of the log-normal. The log likelihood of the best fit decreases steadily with 
increasing Ay because as the constraining power of the data becomes stronger (as reflected in the smaller uncertainties), any deviations 
from the assumed log-normal distribution become apparent. The lower right panel shows Af^ed again, but now color-coded by the number 
of stars in each pixel. Not unexpectedly, at fixed extinction, the uncertainties are higher when the number of stars is lower. 
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5.3.1. Failures of Model Fitting 

In addition to the normal scatter due to random un¬ 
certainties, there are individual pixels where the model 
fitting has obviously failed. These cases are most notice¬ 
able in low extinction regions, where we find occasional 
pixels that have extremely high values of in a re¬ 
gion where there is no other evidence for significant dust. 
There are two major failure modes we have identified in 
these regions. 


Spurious Photometry & the Noise Model — The first of 
these cases corresponds to failures of our adopted noise 
model. If an individual pixel has unusually high contam¬ 
ination from bad photometry, then there will be a larger 
than expected number of spurious detections. When 
these sources are not well-represented by our noise model 
the fitting code will adjust the log-normal parameters 
to increase the likelihood of finding stars far redward of 
the RGB, where such spurious sources typically he (i.e., 
close to the photometric limit). As a result, the best 

fit value of Ay and/or a will be driven to large values 
in an attempt to produce an extremely broad log-normal 
component, increasing the probability of finding very red 
sources. This failure represents a complete breakdown 
of our model. The effect is especially pernicious in re¬ 
gions with low reddening fractions, where one expects 
few stars redward of the unreddened RGB. It is further 
compounded in the outer disk, where the absolute num¬ 
ber of stars per pixel is small. Luckily, however, these 
same failure modes are also marked by much higher than 
average uncertainties in Ay, allowing them to be easily 
identified. An adaptive pixel size would help alleviate 
this issue, at the expense of complicating comparisons 
with other dust tracers by introducing spatial variations 
in the resolution. 

Although the noise-model failures sometimes affect 
only a single isolated pixel, they also sometimes affect 
multiple adjacent pixels (i.e., spanning > 13"). The 
most common reasons for such highly clustered spurious 
sources are (1) improper rejection of sources on diffrac¬ 
tion spikes around bright stars; and (2) subregions of 
background galaxies. To minimize the impact of these 
bad fits in the following analysis, we use unsharp mask¬ 
ing to identify pixels that are wildly different from their 
neighbors, and then replace those pixels with the local 
median of the neighboring pixels. Our rejection thresh¬ 
old is extremely conservative, and therefore some prob¬ 
lem areas clearly remain even in the cleaned maps. All 
subsequent analysis will use these “cleaned” maps. 


WFC3 Calibration Issues — The second failure of the 
models in low extinction regions is due to the WFG3 /IR 
photometric calibration errors discussed in Sec. |4.1[ In¬ 
spection of the right hand plot in Figure shows that 
there are repeating structures in the map m mean RGB 
colo r, ev en after applying the first order correction in 
Sec. |4.1[ These structures follow a square pattern that 
coincides with the locations of the individual WFG3/IR 
frames that make up the survey data. This pattern sug¬ 
gests that some of the width of the RGB is due to stars 
in one part of the detector consistently appearing redder 
or b luer than stars in another area. As discussed in Sec. 
14.11 such an offset could be due to different accuracies 


in the fiat fields used for the FllOW and/or the F160W 
filters, or (more likely) due to systematic differences in 
the accuracy of the spatially-varying model point spread 
functions used for photometry. Attempts to track down 
the exact calibration issue are on-going. We have reduced 
the effect of these structures on the eventual extinction 
map by introducing a fitting parameter Sc to allow very 
small shifts in the color, but the correction is noisy, par¬ 
ticularly when the fraction of unreddened stars is small. 
Thus, while it helps, residual issues remain. 

The pattern is created because our procedure for cre¬ 
ating a model of the unreddened stars identifies stars 
in regions where the RGB is narrow and blue. Regions 
where the RGB is redder and/or broadened by system¬ 
atic errors in the photometry will therefore be best fit 
by models with higher extinction. When these systemat¬ 
ically redder and bluer regions are consistently found in 
certain parts of the WFG3/IR footprint, they imprint a 
pattern of apparently higher than normal extinction on 
the resulting map. However, because these photometric 
issues are at the level < ±0.02 mag nearly everywhere, 
their effects are only noticeable when the extinction is 
small. We can estimate the amplitude of this effect on 
the Ay map by considering the most extreme case, where 
all the stars in the model unreddened GMD would come 
from portions of the WFG3/IR footprint with systemati¬ 
cally blue colors. In this case, unreddened regions that lie 
in the systematically redder parts of the WFG3/IR foot¬ 
print would by redder by < 0.04 mag, corresponding to 
an apparent, spurious extinction of < 0 .3 mag. This 
amplitude is comparable to the typical random errors at 
low extinction (Figures [iq and. We expect these fea¬ 
tures to eventually be eliminatea in the final version of 
the PHAT photometry, which will use spatially-variable 
PSF models in the NIR. 


5.3.2. Possible Systematie Uneertainties 


The uncertainties derived above are tied specifically 
to our chosen model. If this model is incorrect, sys¬ 
tematic uncertainties can be introduced into our anal¬ 
ysis without being reflected in the uncertainties derived 
from the MGMG fit. We have already addressed several 
of the systematics inherent in co nstructing a model of 
the unreddened GMD (Sec. 4.2.5 i.e., failure to identify 
true zero reddening stars, biases from population gradi¬ 
ents and projection effects, foreground reddening, etc), 
but address a few additional points here. 


Variations in the Reddening Law 
We have assumed that there is a single conversion be¬ 
tween reddeniM and extinction that applies throughout 
M31 (see Sec.[^. This assumption may not hold if there 
are local departures from a standard Ry = 3.1 extinction 
law. In such cases we will convert the observed redden¬ 
ing into the wrong extinction. However, the conversion in 
the NIR is relatively insensitive to the value of Ry . Even 
an extreme value of Ry = 5 leads to an inferred extinc¬ 
tion that is only 15% different than in the Ry =3.1 case, 
which is less than our typical precision. We therefore do 
not anticipate this to be a significant source of error in 
the extinction map. 

The Assumption of a Log-Normal Ay Distribution 
Our model assumes that the distribution of extinction 
in each spatial pixel can be approximated as a single log- 
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Fig. 19. — Same as Figure [l^ but for the high extinction star-forming region in Brick 15. Uncertainties I^AyjAy^ lAala^ and Af^ed 
(clockwise from top left) are as a function of median extinction Ay^ for each pixel in Brick 15, color-coded by the parameter indicated in 
the upper right of eac h pa nel. Given the overall high extinction in this region, the associated uncertainties on the parameters are much 
smaller than in Figure [Ts] 


normal. If the true distribution is more complex, then 
the specific values of parameters are harder to interpret 
and their uncertainties may not be correct. 

The log-normal appears to be an adequate representa¬ 
tion for the majority of sightlines in typical Milky Way 
molecular clouds and the diffuse ISM, but the densest 
star-forming clouds are usually found to h ave an addi 


tional powe r-law at high-extinctions (e.g., Kainulainen 
et al. 1 |2009l ; see the more extensive discussion in Sec. 
3|. In practice, deviations from log-normal distributions 
are not a large concern, given that the power-law tail is 
only found in the most actively star forming clouds, and 
contains only a small fraction of the cloud mass when 
present. The highest column densities occupy a very 
small volume of any molecular cloud, and thus random 
sampling of the column density distribution using back¬ 


ground RGB stars will tend to miss any features that 
occupy only a negligible fraction of the pixel area. The 
only notable effect may be a tendency of fits to favor a 
somewhat broader distribution (i.e., larger a), particu¬ 
larly in regions with large median extinctions. Inspec¬ 
tion of Figures [T6| and [TT] suggest that this combination 
of parameter values is not common. 

Another violation of the log-normal assumption may 
occur when the gas distribution is sufficiently complex 
that there are multiple gas layers along the line of sight. 
For example, two molecular clouds that are spatially co¬ 
incident when viewed in projection could have different 
densities and depths along the line of sight, and thus 
would produce a more complex reddening distribution 
than we have assumed. However, even when multiple 
gas overdensities are superimposed along the line of sight. 
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they are unlikely to all contribute equally to the column 
density distribution, such that the total reddening dis¬ 
tribution is likely to be dominated by the subcomponent 
with the highest column density with the largest fraction 
of stars behind it, potentially leaving the log-normal as 
an adequate fit. 

As a potential test of the suitability of the log-normal 
model, we have color-coded the values of the uncertain¬ 
ties by the log of the posterior probability of the best fit 
model in Figures and There is an obvious trend 
such that the log-normal model becomes systematically 
less likely at higher extinction (although only by ^ 1.5 
dex). This behavior is not alarming, because our data 
clearly have less constraining power at low extinctions; 
when the ability to constrain parameters is low, mak¬ 
ing it impossible to rule against a given adopted model. 
Any model is therefore acceptable in terms of the poste¬ 
rior probability of the best fit. While we could engage in 
more exhaustive testing of whether a log-normal is the 
optimal functional form, we proceed with an understand¬ 
ing that the log-normal parameters are best interpreted 
as a good, but approximate, first-order characterization 
of the distribution of column densities within the pixel. 

Biases from Reddened Young Stars 


Our model assumes that all of the stars that we fit in 
the CMD are reddened RGB stars. However, in regions 
of significant reddening, some of the younger stars that 
are blueward of our analysis region can be reddened into 
the portion of the CMD that we are analyzing. These 
younger st ars are typically found in high extinction re - 
gions (e.g.,|Harris et al.|1997t|Zaritsky et al.|2002||2004|), 
making this mechanism a possible source of bias. 

Unfortunately, the amplitude of the bias is difficult 
to quantify because the young stellar populations have 
much more dramatic spatial variations (Lewis et al. 2015) 
than the older RGB stars, which have been well-mixed 
over many dynamical times 250Myrs in M31’s main 
star-forming ring). Thus, the exact degree to which 
highly reddenened young stars contaminate the RGB 
analysis region will depend on exactly how numerous the 
young stars are relative to RGB stars. The degree of 
contamination will also depend on the ages of the young 
stars, since the core Helium-burning sequences (roughly 
25-400 Myr) are redder than the main sequence, and thus 
can more easily be reddened into our analysis regions. 
Finally, the degree of contamination will also depend on 
the depth of the observations. At brighter magnitudes, 
core Helium-burning stars are well-separated from the 
RGB, but at fainter magnitudes, they merge into the 
red clump where their colors are indistinguishable from 
the main population of RGB stars. In this latter case, 
reddening of these stars would produce no bias, because 
they would be included in the unreddened model. 

Although the exact amount of bias is difficult to com¬ 
pute, we can make a first order estimate by calculating 
the colors of stars at different ages, and their relative 
numbers compared to the RGB. We adopt photometric 
errors comparable to those in PH AT, and assume a con¬ 
stant star formation rate. Overall, this choice will tend 
to overestimate the contribution of young stars, given 
that M31’s st ar formation rate has declined since a burst 
2- 3 Gyr ago (iWilliams et al.||2Q15|). However, Figure 13 
of Lewis et ar~( ^bl5| ) shows that there are some star- 
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Fig. 20.— The distribution of colors for stars in 4 different age 
ranges (100-500 Myr, 500Myr-lGyr, 1-2 Gyr, and >2 Gyr; black, 
blue, magenta, and red, respectively), assuming a modest degree 
of metallicity evolution ([Fe/H]= 0 for the two most recent bins, 
-0.25 for the 1-2 Gyr bin, and —0.75 in the oldest age bin), for a 
simulated GMD assuming a constant star formation history and 
typical PHAT photometric errors. Each panel shows the distribu¬ 
tion of colors within a 1 mag interval of F160W magnitude. To 
assess the possible contribution of reddened young stars to the 
RGB-dominated portion of the GMD, all colors are given relative 
to the peak of the stars in the oldest age range; the arrow in the 
bottom left panel shows the reddenening expected for Ay = 1 mag 
of extinction. The numbers listed in each panel give the number of 
stars in each histogram; percentages are given relative to the num¬ 
ber of stars in the oldest peak. Stars older than ~ 500 Myr have 
NIR colors that are indistinguishable from the main RGB peak. 
Stars in the youngest bin can have bluer colors at young ages, but 
are not numerous and require large reddening before they could be 
mistaken as RGB stars with modest reddening. 


forming regions that have recent star formation rates 
(averaged over 100 Myr and ^ 0.4 Gyr timescales) that 
are of the order of the past average star formation rate 
or higher. In these regions, the potential bias will be 
somewhat higher than we estimate here. 

In Figure we plot the resulting distribution of col¬ 
ors for 4 dinS*ent age ranges (100-500 Myr, 500 Myr- 
1 Gyr, 1-2 Gyr, and > 2 Gyr), assuming a modest de¬ 
gree of metallicity evolution ([Fe/H] = 0 for the two 
most recent bins, -0.25 for the 1-2 Gyr bin, and —0.75 
in the oldest ag e bin); the firs t age r anges encompass the 
ages probed in Lewis et al. (2015)’s spatially resolved 
maps of M31’s recent star formation history. We do 
not include results for ages younger than 100 Myr, since 
these stars are too few and too blue to make any signif¬ 
icant contribution. We plot the distributions relative to 
the color of the oldest stellar populations (i.e., relative 
to the peak of the RGB), for four different magnitude 
ranges. Each plot lists the numbers of stars redder than 
FI low — F160IU = 0.2, both as an absolute number 
and as a percentage of the stars in the oldest RGB peak; 
these numbers give some sense of the relative size of the 
young population that might potentially contaminate the 
RGB when reddened. 
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Figure suggests that the impact of contamination 
from young reddened stars is likely to be small. Stellar 
populations older than ^0.5 Gyr have colors that are 
indistinguishable from the RGB. While their presence or 
absence may affect the total number of stars in a given 
pixel, our model fits only relative distribution of stars 
on the GMD, and is thus insensitive to their presence. 
Stars in the 100-500 Myr age range will likewise have no 
impact at faint magnitudes (F160W> 22) where their 
fractional contribution is small (< 5% of the number of 
RGB stars), and their colors are either indistinguishable 
from the RGB (22 < F160W < 23), or too blue to make 
any contribution except at the very highest extinctions 
(23 < F160VF < 24). 

At brighter magnitudes (22 <F160W), the 100- 

500 Myr stars are proportionally more important. They 
contain up to 20-25% of the number of stars found in 
the RGB (assuming a constant star formation rate), and 
could be reddenened to the color of the RGB peak with 
1-1.5 mag of visual extinction. If the young stars ex¬ 
perienced this amount of extinction, they would merge 
with the peak of unreddened RGB stars. The model fits 
would then assume that the fraction of reddened stars 
was smaller than it actually was. To affect the measure¬ 
ment of Ay, however, the stars would need to experience 
even more reddening, such that they moved redwards of 
the RGB peak. In this case, the presence of the stars 
would tend to bias the measurement of Ay to lower val¬ 
ues. However, given that (1) RGB stars would still out¬ 
number the young stars 3:1 at bright magnitudes; (2) 
that the fainter stars where the bias is low are factors of 
4-6 times more numerous; and (3) that the inclusion of 
the noise model can potentially “absorb” some fraction 
of stars that appear abberant compared to the pure RGB 
fit, it is not clear by how much the best fit model would 
be biased in practice. We therefore conclude that biases 
due to reddening from young <0.5 Gyr old stars is likely 
to be negligible in most cases, and would only potentially 
be detectable at the highest extinctions {Ay > 2), at 
anomalously high SFRs in the 100-500 Myr age range, 
and where the data was shallowest (i.e., in the inner 
disk). 

Loss of Stars at High Ay 

Another possible source of bias would be if highly red¬ 
dened stars are systematically missing from our data in 
high extinction regions, due to their suffering enough ex¬ 
tinction that they fall below our detection limit. Given 
the range of extinction found in our analysis, rarely is the 
extinction high enough to completely remove large num¬ 
bers of moderately bright RGB stars (e.g.. Figure [^. 
This is not to say that no individual stars have been r^- 
dened out of the sample. Intrinsically faint RGB stars are 
certainly “missing”, and undoubtedly occasional bright 
RGB stars may be extinguished as well if they happen to 
fall behind the very smallest, densest core of a molecular 
cloud. However, it appears highly unlikely that either 
effect is not being accounted for in the model fitting, as 
long as the mode of the log-normal distribution is de¬ 
tectable for some of the RGB. 


Summary of Systematic Biases 

All of the poten tial systematic biases identified above 
and in Sec. 4.2.5 are modest (Ay < 0.3 mag) or con¬ 
fined to very rare regions. Systematic errors are therefore 


likely to only be noticeable at low extinctions. However, 
the random errors at these low extinctions are large, so 
that the systematic errors may not actually dominate the 
error budget for any individual pixel. However, the ef¬ 
fects of systematics can be noticeable when looking at 
the cumulative behavior of the extinction over large, low 
extinction regions of the galaxy, which essentially aver¬ 
ages over many pixels and thus reduces the impact of 
random unce rtain ties. We look for evidence of these ef¬ 
fects in Sec. |6.4| below, after presenting our final dust 
maps and comparing them with other measurements of 
dust extinction. 


6. RESULTS 


6.1. A Galaxy Wide Map of M3Hs Cool ISM 
In Figur e!^ we present the full maps of M3Fs median 
extinction Ay. Figures 


22 


24 


show a series of zoomed 
in regions of the same maps^Rom the inner to the outer 
disk (i.e., moving from bottom to top in Figures [21]). All 
images show the full interleaved map (12.5 pc pixel sam¬ 
pling with 25 pc resolution). We present only the cleaned 
maps; some remaining spurious features due to high num¬ 
bers of false detections near bright stars are still visible, 
particularly in the outer dislp^ There are also very low 
level “grid” structures due to the ±0.01 mag structured 
systematics in the p hotom etr y acro ss the WFC3/IR chip 
(as discussed in Sec. |4.1| an d |5.3.2| ). 

The maps in Figures jzTj—124| reveal extremely rich de¬ 
tails. Assuming that the extinction is proportional to the 
column density of cold gas, the maps of Ay offer one of 
the highest resolution views of the ISM of a spiral disk, 
outside our own Milky Way. The structure of the ex¬ 
tinction shows networks of filaments (Figures fc 24), 
holes (Figure]^, ripples/feathers (Figure]^ seen as re¬ 
peating hook^ structures on the inside edge of major 
star-forming ring), and spurs (Figure]^ linear features 
extending radially from the major starTorming rings and 
spiral arms). These latter two features are thoug ht to be 
direct results of hydrodynamic instabilities (e.g., Kim fc 


Ostriker|2002l[C 

20061 IKim & Oi 


Ghakrabarti et al. 


|2003[ IShetty 

)bs & Bonnel 


^ 0 striker 


20061 |Kim fc 6striker||2006[ |Dobbs ^ 


ll||2006l|L ( 


ee 


& Shu|2012), making the maps a superb resource for un- 

derstanding these common features of s piral disks (e.g., 
Elmegreen 1980 La Vigne et al.||2006). For example, 
the regular pattern of spurs evident m the inside edge 
of M3Fs star-forming ring (Figure [2^ has a particularly 
small opening angle, and quickly bends around to form 
almost a complete inner ring. In contrast, simulations of 
spur formation typically find much more open patterns 
of spurs than are seen here. 

Inspection of Figure [^ shows that typical extinctions 
are of order Ay ^ 1, with very little surface area hav¬ 
ing Ay >3. If we limit the map to only those re¬ 
gions with well-measured values of Ay (fractional er¬ 
rors A Ay/Ay <0.3), the median of Ay is 1.0 mag, with 
only 10% of the pixels having Ay >1.8, and 1% having 


There are also som e slight problems visibl e in Field 8 of Brick 
22 (for nomenclature see|Dalcanton et al.|2012|) that can be traced 
to unresolved photometric issues. i'Jianktully, these are in very low 
extinction regions and thus the scientific impact of these problems 
is low. There is also an issue in Field 12 of Brick 3, wherein the 
photometry is systematically too red. 
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Fig. 21.— Map of the median extinction Ay for all analyzed regions. Some very low-level “grid” patterns are visible in the map, due to 
position-dependent calibration issues with the WFC3/IR chip (Sections |4.1|&: [5.3|>. 


Ay > 2.8. For comparison, the Milky Way sample o f 23 
molecular clouds studied by Kainulainen et al. (2009) has 


a median of Ay = 1.2, and only 2 clouds have Ay >3. 
Thus, although one cannot draw fi rm statistical con¬ 
clusion s based on the heterogeneous [Kainulainen et ^ 
(2009) sample, the median extinctions in M31 are not 
obviously atypical for what is seen locally in the Milky 
Way. 

One way to assess the veracity of the extinction maps 
is to compare to other tracers of extinction. A straight- 
foward test is to compare the dust morphology to that 
seen in optical imaging in blue filters. At these wave¬ 
lengths, dusty structures can appear as dark features 
against the smooth background of unresolved stars, al¬ 
beit with confusion from foreground stars. Optical im¬ 


ages also have the highest spatial resolution of various 
extinction tracers. 

In the lower rig ht p anels of Figures [25]—[Tf| we show 
blue optical imagej^ of the same reg ions as the extinc¬ 


tion maps in Figures 


22 


(reproduced in the 


23 and 24 

upper left of Figures 25 27\ for reference). The extinc¬ 
tion features seen in our maps rival the resolution of the 
equivalent features seen in the optical images. The mor¬ 
phology of the features are in excellent agreement as well. 

As a further morphological comparison, the lower left 
panels of Figures [25]— sh ow the publicly availabl e 
Westerbork maps oi HiRom Brinks & Bajaja (1986). 


Courtesy of the superb astrophotographer Robert Gendler: 
http://WWW.robgendlerastropics.com/ 
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Fig. 22. — Map of the median extinction Ay for the inner regions containing Bricks 2, 4, 6, and 8 along the lOkpc star-forming ring 
(bottom to top), and the major axis field Brick 5, which samples a smaller, less intense star-forming ring. Note the regular ripples or spurs 
along the inside edge of the star-forming ring. 


Although these maps have significantly lower resolution 
(45") than the extinction maps, the broad morphology 
is again in good agreement. 

While the comparison between the optical image and 
the extinction map show many of the same detailed struc¬ 
tures, the new extinction maps offer a far less ambigu¬ 
ous view. First, they lack the complicated foreground of 
young stars that confuses interpretation of extinction fea¬ 
tures in broad-band images. But more importantly, the 
new maps measure the total extinction of the dust col¬ 
umn, regardless of where it is with respect to the stars. In 
contrast, there are extinction features in the broad-band 
optical image which appear extremely strong, when in 
fact the strength results only from the dust layer being 


closer to the foreground, allowing it to block more of light 
due to geometry alone. A clear example of this effect can 
be seen in the optical image in Figurej^ where the struc¬ 
ture at (10.92°, 41.42°) appears to be much higher col¬ 
umn density than the neighboring structure at (10.99°, 

41.45°). Comparing to the adjacent map of Ay, how¬ 
ever, reveals that the latter clump is actually the higher 
extinction region. Instead, the high apparent extinction 
is due primarily to the much higher reddening fraction 
in that region. This example points to the peril of us¬ 
ing color and/or unsharp masking to infer the extinction 
from maps of the attenuation (e.g., Regan et al. 2011), 
particularly in inclined galaxies with large variations m 

fred' 
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Fig. 23.— Map of the median extinction Ay for: (1) the major axis of the the lOkpc star-forming ring (Bricks 12, 14, and 16 along the 
left, and Bricks 15 and 17 along the major axis, from bottom to top); (2) the major axis field Brick 9, which samples the second of the 
2 inner star-forming rings; and (3) the star of the outermost star-forming ring (Brick 18, in the upper left). The ripples/spurs noted in 
Figure [^continue into Brick 12 on this map. 

6.2. Comparison with Emission-Based Maps of Dust 

Mass 

Thanks to revolutionary new mid- and far-infrared fa¬ 
cilities ( e.g., Spitzer, Herschel), and sophisticated mod- 
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Desert et al.|l990l Draine & Lij2UUY||Uompiegne 
now become standard to use long- 
wavelength spectral energy distributions (SEDs) to in¬ 
fer the mass, composition, and temperature of the dust, 


a function of position, which, for a given dust model, can 
be converte d directly into pred ictions for the extinction. 


and the illuminating radiation field 
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Recently, Draine et al. (2014) published a state-of-the- 
art analysis of lV131's dust, using emission maps from 
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Fig. 24. — Map of the median extinction Ay in the outer disk, including Bricks 17, 19, 21, and 23 along the major axis, and Bricks 18 and 
22 on the left of the image. Two of the fields in Brick 22 (near a = 11.86°, (5 = 42.14°) have somewhat higher-than-expected photometry 
errors, which leads to a large fraction of erroneous pixels in the dust map. 


r^BOOjum, using observations in 13 different bandpasses. 
The resulting maps of the dust column density were gen¬ 
erated at two resolutions: 24.9" and 39", corresponding 
to the resolution of the SPIRE 350/im and the MIPS 
160/im observations, respectively, and then sampled with 
either 10" or 16" pixels. A high (6") resolution emission 
map has also been published by Montalto et al. (2009), 
but the level of contamination from stellar sources makes 

it difficult to compare to the work here. _ 

The emission-based dust map of Draine et al. (2014) 
can be used to make a direct prediction tor the extinction. 


Draine et al. (2014) derive the extinction-to-dust mass 
surface density ratio AyjYjMd = 7.394 for their specific 
model of the dust composition, where T^Md is the dust 


mass surface density in Mq/ pc^. 

In the upper right of Figures [25] —[27| we plot the pre¬ 
dicted distribution of^^ deri^^ Rom the 24.9" reso¬ 
lution Draine et al. (2014) dust map (i.e.. Ay, 

The niorphological agreement between the Ay^ 


,ermsston 


,ermsston 


map and the CMD-based Ay map is superb, particu¬ 
larly in regions where Ay > 1 mag. Even very small 
features are clearly reproduced in both of the maps, in 
spite of their completely independent derivations. It is 
hard not to see the level of agreement as a triumph for 
both techniques. 

That said, the maps do present some differences. The 
CMD extinction maps have nearly 4 times the spatial res- 
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Fig. 25.— Comparison between the new M31 extinction map ( Ay ; upper left), the extinction inferred from the emission-based Drai ne e t 
al dust map (upper right), an Hi map (lower left), and an optical image (lower right). These maps reproduce the area shown in Figure [2^ 


olution (93 pc vs 25 pc), which naturally produces clearer 
views of the smallest features. However, the emission- 
based maps will be sensitive to any dust in a given pixel, 
whereas the CMD extinction-based maps can miss dust 
that is locked into small dense structures with negligible 
filling factors. The emission-based maps are also possi¬ 
bly more accurate at low extinctions, where the CMD- 
mapping technique is limited by our ability to: (1) detect 
broadening over the intrinsic width of the RGB; and (2) 
generate accurate “zero-reddening” models of the RGB. 
On the other hand, the emission-based maps face their 
own difficulties at low extinctions, where the uncertainty 
in wide-field background subtraction and calibration can 


easily produce large-scale systematic offsets in the pho¬ 
tometry. Emission maps may also have difficulty disen¬ 
tangling emission from unres olved dust componen ts with 
different temperatures (e.g., Galliano et al. |2011 ). 

The most glaring discrepancy between the two tech¬ 
niques is the offset in the normalizations. Examina¬ 
tion of the color scales in Eigures [25]— shows that the 
emission-based maps predict extinctions that are more 
than twice those that are observecP^ This difference is 
present even in low extinction regions, suggesting that it 


Note that this offset is not due the RGB only sampling half of 
the dust layer, because we explicitly fit for the fraction of reddened 
stars (eqn.T^ when modeling the CMD. 
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Fig. 26.— Comparison between the new M31 extinction map ( Ay ; upper left), the extinction inferred from the emission-based Drai ne e t 
al dust map (upper right), an Hi map (lower left), and an optical image (lower right). These maps reproduce the area shown in Figure [2^ 


cannot be due to simply missing very dense dust clouds 
with very small filling factors, which would only affect 
the highest column density regions. Given the simplicity 
of the CMD-based extinction technique compared to the 
emission based technique (which relies on assumptions 
about the radiation field as well as calibrations w ithin 
the Milky Way; e.g., Weingartner & Draine 2001), it is 
far more likely that tnis global oftset represents an issue 
with the current calibration of the dust emission models. 

We make a more quantitative assessment of this off¬ 
set by comparing a resolution-matched version of our 
CMD-based extinction map to the emission-based map. 
We use the mean extinction (Ay) for all comparisions. 


rather than the median Ay; these two quantities can 
be significantly different for a log-normal distribution, 
and the mean is the better estimate of the total dust 
column density and thus is a fairer test of the quan¬ 
tity tracked by the emission-based map. We extract 
the posterior probability distribution function of (Ay) 
directly from the MCMC samples, rather than deriv¬ 
ing it from the combination of the best fit Ay and 
a (i.e., (Ay) = exp (cr^/2); Eqn. [^; this choice 

avoids any biases due to correlations between Ay and 
a. We then degrade the resolution of the {Ay) map 
from FWHM = 6.645" to FWHM = 24.9" by smoothing 
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Fig. 27.— Comparison between the new M31 extinction map (Ay; upper left), the extinction inferred from the emission-based Drai ne e t 
al dust map (upper right), an Hi map (lower left), and an optical image (lower right). These maps reproduce the area shown in Figure [24] 


with a Gaussian ke rnel of st andard deviation of 10.19" 
(= \/24.9^ — 6.645^/2V21n2); this degrades the resolu¬ 
tion from 25 pc to 94 pc. We use this resolution-matched 
map of (Ay) in all subse quent comparisons with results 

from Draine et al. (2014). _ 

The top left panel in Fig ure shows the p ixel-by- 
pixel correlation between the Draine et al. (2014) extinc¬ 
tion maps and the mean extinction {Ay) of the reddened 
component, derived from the NIR CMD. The pixels are 
highly correlated, as expected from the excellent mor¬ 
phological agreement seen in Figures [25|—[27] However, 
it is clear from t he top left of Figure [28 [that the ex¬ 
tinctions from the Draine et al. (2014) maps are higher 


by a factor of >2. There is also a small tail of points 
where there is little dust emission, but the CMD-based 
extinction is higher. These points are all due to regions 
where the photometry is slightly biased to redder colors 
near the chip edge in the very outer disk, where there 
are other photometry errors (Fields 2 and 8 of Brick 22), 
or where crowding is high and the RGB is intrinsically 
broad in the very inner disk (Fields 3 and 4 of Brick 3). 

6.2.1. Correcting the Emission- and CMD-hased Extinctions 

We quantify the differences between the emission- 
based and CMD-based extinction maps by calculating 
a transformation that would bring the maps into agree- 
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Fig. 28.— (Top Row) Pixel-by-pixel comparisons between tlie|Draine et al.|(|2014|) emission-based Ay map and the resolution-matched 
map of mean extinction (Ay) derived from the NIR CMD. The original comparison (left) shows strong correlations, as expected from the 
morphological agreement. However, the steep slope and the vertical offset from the origin indicate issues with the overall sc aling between 
the two methods, and a bias in the CMD-based extinction maps at low dust column. After correcting the|Draine et al.|(|2014|) maps (right) 
by an overall scaling factor of R = 2.53, and correcting the CMD-based extinction by an additive bias faHoF'oFS'^^Qn^T^he agreement 
between the two maps is much improved. The scatter is still significant, however. (Bottom Left) The pixel-by-pixel ratio between the two 
maps and the resolution-matched map of mean extinction (Ay) derived from the NIR CMD. The red line shows the predicted ratio using 
the scale+bias correction used in the center panel. The results are similar if the median extinction of the reddened component is used 
instead of the mea n, although the scale factors are slightly different. (Bottom Right) The distribution of the values of the scale factor 

correction R to the 


Draine et al. 


( 20i 4 map and the bias correction b to the CMD-based map, derived where Ay > 0.5. 


ment. We assume that the Draine et al. (2014) extinc¬ 
tions are too high, and must be divided by a lactor R. 
We include an additional correction factor b for possible 
systematic bia ses in the CMD-based extinctions (see Sec. 


4.2.5 & 5.3.2). We then estimate appropriate values for 
R and b by calculating on n grid of R and 6, where 
is calculated using the observed ratio Ay^ernission/{^v) 


compared to the prediction for the model as a function of 


{Ay) (i.e., Ay^emission/(Ay) = R{l^b/{Ay))] lowcr left 
of Figure!^. When calculating we use an empirical 
estimate oTthe uncertainty in the measured ratio as a 
function of (Ay), since we lack a realistic noise model 
for the ratio R. We limit the calculation to only points 
with Ay > 0.5 mag where measurements for both dust 
models are expected to be more robust. 

The resulting distribution is shown in the lower 
right panel of Figure The best-fit model parame¬ 
ters are R = 2.53 and b = —0.04 mag when using the 
mean CMD-based extinction {Ay). If the CMD-based 
extinctions are assumed to be unbiased (5 = 0) then the 


scaling for the Draine et al. (2014) extinctions shifts to 
slightly lower values of R. 

In the upper right panel of Figure we show 
the pixel-by-pixel correlation after applying these cor¬ 
rections, such that emission ^ ^F,emission/ and 
{Ay) -A {Ay) + b. The prediction for the ratio using 
the same model is also shown as the red line on the 
lower left panel. The overall distribution now has the 
desired slope of 1 (red line) and is consistent with pass¬ 
ing through (0,0). 

The fits shown in Figure 28 suggest that (1) there is 
a negligible global bias in tne CMD-derived extinction 
values, corresponding to a color shift of ^ 0.005 magni¬ 
tudes in FllOVF—F1601F, and (2) there is a large scaling 
difference between the emission- and CMD-based extinc¬ 
tions. We now discuss the second of these conclusions in 
mor e de tail, and defer the less pressing issue of biases to 
Sec. 16.41 
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6.2.2. Are the Draine et al. 1(2014]) Extinctions Really Too 

- 

Figure strongly suggests that there is a factor of 
^ 2.5 difference between the extinction inferred from 
CMD fitting and that inferred from modeling the IR 


emission. The distribution of values in Figure |28 
shows there is no reasonable value of the bias for which 
the ratio R becomes 1. Thus, we are in the uncomfort¬ 
able position of accepting a large systematic problem in 
either the CMD-based extinction or in widely used mod¬ 
els of dust emission. 

We argue that the factor of 2.5 scaling difference is 
more likely t o have resulted from a problem with the 


Draine et al. (2014) emission-based maps. Simply put, 
there is no obvious way that the extinction inferred from 
the CMD can have mistaken the degree of reddening by 
such a large factor. We demonstrate this in Figure 
which shows the CMD of the same region plotted in Fig¬ 
ure [m The underlying density map shows the CMD 
expectation for the true best fit model (upper left) and 
for a model where the true reddening is ^2.5 times higher 
than what we derived (shown in the lower row as an ad¬ 
justment to the median extinction Ay or to the width of 
the distribution, on the left and right, respectively). A 
higher extinction is completely at odds with the actual 
observations. 

The only ways out of this conclusion are deeply un¬ 
satisfying. For the CMD-based extinctions to be artifi¬ 
cially low at essentially all pixels, we would have had to 
(1) mistaken the true location of the unreddened CMD 
by many tenths of a magnitude, or (2) invoked a red¬ 
dening law with a drastically different relationship be¬ 
tween attenuation and reddening in the NIR. We can 
rule out the first escape route because we have many re¬ 
gions where the vast majority of stars are unreddened 
(i.e., small fred)] in these regions, the NIR RGB is es¬ 
sentially exactly where our model places it, and is fully 
consistent with theoretical isochrones. Moreover, such a 
color offset would manifest as a constant bias factor, not 
a difference in scaling. The second solution is equally 
difficult to accept. While r eddening laws are k nown to 
vary in the ultraviolet (e.g., Gordon et al.|2003), there is 
no evidence for significant variations in the NIR at the 
required amplitude. 

In contrast to the difficulty explaining a large scale er¬ 
ror in the CMD-based extinctions, there are many pos¬ 
sible ways to explain scale errors in the emission-based 
extinction maps. Converting the spectral energy distri¬ 
bution of the dust emission to a prediction for the ex¬ 
tinction requires calculating the dust mass in the con¬ 
text of a specific model for the physical properties of the 
dust (composition, grain size distr ibution, etc) and for 


the starlight radiation held. The Draine e t al. (12014 


dust map is derived using the Draine Li "(2007) mo 
els, which are calibrated to niatch the eni ission of dif- 
fuse, dusty regi ons in the Milky Way (e.g., Weingartner 


& Draine 2001), with additional adjustments in the vol¬ 
ume of gram material to give better consistency with 
deplet ion measurements o f the dust mass per hydrogen 
atom (Draine et al. |2014). Given the complexity of the 
models and the calibration, it is possible that these lo¬ 
cally calibrated models do not apply in a different envi¬ 
ronment. Indeed, Draine et al. (2014) specifically raises 


the possibility of their systematically overestimating the 
dust mass when discussing differences with Planck emis¬ 
sion, although they argue that systematic errors initially 
appear smaller than a factor of ^1.3. 

Additional complications arise from the degeneracy in¬ 
herent in fitting dust models to spectral energy distri¬ 
butions, and from the inability to fully re solve complex, 
multi -temperature dust distributions (e.g., Galliano et al. 


2011). Dust emission is highly non-linear with temper 


ature, and thus the ensemble emission from a collection 
of dust clouds with varying temperatures does not ac¬ 
curately reflect the mean mass-weighted temperature. 
There is also evidence t hat cold dust may b e much more 
emissive than expected (jParadis et al.|2009|) . If the com¬ 
ponent of high-emissivity dust is pervasive in M31 as well 
as the MW ^ then it may help to explain why the Draine 


et al. (2014) map infers more extinction than is actually 
seen. However, this effect seems restricted to the molec¬ 
ular phase, which is not dominant in M31. 

Gomparable mism atches between extinction and the 


Draine & Li (2007) models have been s een b efore in 
the literature. Both Alton et al. (2004) and Dasyra 


et al. (2005) found that multiwavelength radiative trans- 
fer niodeling of edge-on galaxies underpredicts the ob¬ 


served sub- mm emission at 850/im when using the |Draine 


& Li (2007) models. In other words, to be consistent with 
the observed amount of extinction, either there must be 
another component of dust not accounted for in their ra¬ 
diative transfer model, or the dust must produce a factor 
of ^3-4 X gre ater emission at long wavelengths than pre¬ 
dicted by the Draine & Li ([2007) model. The inverse of 
this latter possibility — i.e., that a given amount emis¬ 
sion can be produced by 1/3 the amount of dust - would 
be consisten t with our resu l ts, aft er including the model 
revisions in Draine et al. (2014). Recent simulations. 


however, suggest that the mismatch may instead be due 
to asymmetries and inhomogeneities in the disk that af- 
fect the outcome of radiative transfer calculations (e.g. 
Saftly et al.||2015 ) 


Even mor e clearly, recent work by [Plank Gollaboration 


et al. (2014) showed a nearly identical discrepancy to the 
one we identify here. They measured the optical extinc¬ 
tion observed towards thousands of QSOs observed with 
the Sloan Digital Sky Survey and compared it to the 
extinction predicted by fitting ([Draine & Li[[2007[) mod¬ 
els to all-sky data from Pla nck, IKAE, anR ITlb A'. This 
comparison found that the Draine & Li (2007) extinc¬ 
tion estimates were too high by a factor of rJ 2 — 2.4. 
They found a comparable discrepancy when looking at 
molecular clouds whose NIR extinctions were mapped 
using 2MASS by [Sch neider et al.[ ([2011[). We compare 
our results to those m Plank Collaboration et al. (2014) 
in more d etail in Sec. [6.5 [ below. 

Finally, Lombardi et al. (2014) find a factor of ^2 dif- 


ference between obse rvations in Orion and the [Weingart-[ 
ner & Draine ([2001 ) predictions for the ratio of 2.2/im 
extinction to the 850/im opacity. They arg ue that this 
discrep ancy is reduced with recent models by [Ormel et al.[ 
( 2011|. Interestingly, however, the ratio measured by 


Plank Gollaborati on et ah ([2014[) is actual l y in m uch bet- 
ter agreement with Weingartner & Draine (2001), in spite 


of showing a siginihcant difference between the measured 
an d the p redicted extinctions from the similar Draine & 
Li] ([2007) models. 
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Fig. 29. — Binned CMD showing the reddened model derived for stars found in a 25 pc pixel in Brick 15 (as in the right panel of 
Figure pT] ). The upper row shows the actual best fit model derived from CMD fitting (left) and the joint probability distribution f or (Ay) 
and cr (right j. The lower row shows example models corresponding to a factor of 2.53x higher mean extinction, as inferred from the|Draine| 
|et al.|(|2014| dust emission model, made by rescaling the median extinction (left) or the width of the reddening distribution (right). Clearly, 
tlie data fsnown as blue stars) are not consistent with the higher mean extinction characteristic of the emission-based dust model. Note: 
slight differences in the best fit parameter values given on each plot are due to difference instances of fitting the same pixel. 


In summary, it appears that the most widely used 


Draine & Id] (|2007|) dust models overestimate the mass 
ot dust needed to produce a given amount of emission 
in M31, even after revisi ng the dust masses d ownwards 
using the modifications in Draine et al. (|2014). This dis¬ 
crepancy is particularly troubling, given the ever growing 
importance of dust modeling in the era of flagship mid- 
and far-IR observatories. Further application of CMD- 
based extinction measurements outside the Milky Way 
will be invaluable in improving existing dust models. 

6.3. Small-Scale Variations in Emission- and 
CMD-based Extinctions: The Impact of the 
Interstellar Radiation Eield 

In Figure [3Q|we plot the ratio b etween the emission- 
based Ay from Draine et al. (2014) and the CMD-based 


mean extinction (Ay), after smoothing the latter to 
match the resolution of the former. The left panel plots 
the raw ratio, and the right pane l plots the ratio after 
correcting the Draine et al. (2014) extinctions by a fac¬ 
tor of 77 = 2.53 and the CMD-based mean extinctions 
by an offset b = —0.04 mag. The dark contour shows 
the locus within which {Ay) > 1 in the smoothed map; 
CMD-based extinctions outside this locus will potentially 
be subject to larger systematic uncertainties. The red¬ 
der regions of the right han d panel indicate regio ns where 
the mismatch between the Draine et al.| (2014) and the 
CMD-based extinctions are larger. In the bluer regions 
the agreement is better, but only approaches good agree¬ 
ment emission/ {^v) ^ 1 in the left panel) at the very 
lowest extinctions, where both methods agree there is 
very little dust. 

































Dust Mapping in M31 


39 


Figure 30 shows an amount of scatter consistent with 
that seen in the pixel-by-pixel correlations shown in Fig- 
ure|^ However, it also shows that the deviations from 
the mean ra tio a re often spatially coherent. In the follow¬ 
ing section (6.4) we assess the possible contribution from 


systematic errors in the CMD-based extinctions, but be¬ 
fore doing so we discuss how such behavior could also be 
produced by emission-based extinction measurements. 

Non-linearities in Smoothing — To produce Figure |3Q| we 
have smoothed the CMD-based dust map. However, this 
smoothing is weighted by the column-density of the dust. 
In contrast, the smoothing that takes place during emis¬ 
sion observations is flux-weighted. The emission from 
dust is a strong power of dust temperature, and thus the 
resolution of a mid- or far-IR telescope will effectively 
smooth the dust map non-linearly, with more weighting 
given to regions with higher dust temperatures. This can 
lead to em ission maps missing cold dust (e.g., Galliano 
et al.|2011), and producing systematic mismatches with 
the observed extinction when there are strong, sharp gra¬ 
dients in the dust temperature. This latter situation is 
likely to be found on the edges of star-forming regions, 
such as just beyond the boundaries of the star-forming 
ring (see Figure 9 of Draine et al.||2Q14~ ). Such an effect 
could be responsible for the more significant deviations 
(redder colors) frequently seen just outside the {Ay) = 1 
contour. 


Variations in Dust Composition — Models of dust emis¬ 
sion typically allow the dust composition, and in par¬ 
ticular the poly yomatic hydrocarbo n (PAH) content, to 
vary spatially. Draine et al. (2014) fit for PAH abun¬ 
dance Qpah and find strong spatial variations in and out 
of M31’s star-forming rings, as well as large gradients 
within the galaxy. Some of the abundance variations are 
lessened when variations in the interstellar radiation field 
are also included, particularly in the inner galaxy, but the 
change in PAH abundance across major star-forming fea¬ 
tures appears to persist in the more structured northern 
half of the galaxy. Interestingly, the largest red region 
outside the (Ay) = 1 contour is also a region where the 
PAH abundance appears anomalousl y high for such a low 
column density region (Figure 11 of Draine et al.||2014). 
It may be possible that coherent variations in the inferred 
dust composition may be propagating into coherent vari¬ 
ations in the dust extinction. It would not be surprising 
if M31’s dust composition varied systematically between 
the diffuse gas in the low column density regions and the 
denser gas in high column density regions. It is also pos¬ 
sible that some of the bias factor b that we have ascribed 
to the CMD-fitting method may instead be a systematic 
in how the dust composition changes in the more diffuse 
ISM. We re turn to the issue of possible PAH correlations 
in Sec. 


The Speetrum of the Interstellar Radiati on Field — One of 
the most interesting features in Figure [3Q| is the presence 
of two large regions at high extinction where the ratio 
Ay ^emission I {Ay) is much lowcr than typical - one in the 
lOkpc star-forming ring just west of the major axis, and 
one in the outermost spiral arm along the major axis. 
These two regions are the locations of two well-known 
recent star f ormation events. Th ey can clearly be seen in 
Figure 13 of Lewis et al. (2015), which plots the “Scalo 


6” parameter, defined as the current star formation rate 
{SFR) divided by the past average SFR. This parame¬ 
ter is closely associated with the hardness of the under¬ 
lying spectrum. As the Scalo b parameter increases, so 
does the fraction of young stars, which will increase the 
contribution from hot stars to the underlying interstel¬ 
lar radiation field, potentially changing the efficiency of 
dust heating for a given bolometric energy density in the 
interstellar radiation field. 

To explore a possible empirical correlation between the 
spectrum of the underlying stars and the spatial varia- 
tion of Av,emissionl{Av), in Figure 


31 


we plot the ob¬ 
served ratio of dust extinctions as a function of the Scalo 
b parameter (= SFRiqq/{SFR)) where the “current” 
star formation rate SFRiqq is integrated over the past 
100 Myr using the spat ially-resolved recent star forma¬ 
tion history derived in Lewis et al. (2015), derived on 
100 pc bins. We also plot a rolling average constructed 
by taking the mean of 200 sequential points ranked by 
the Scalo b parameter. We restrict this comparison to 
regions where {Ay) > 1, which should better sample the 
denser neutral and molecular ISM. 

Figure reveals that there is a systematic drop in 
Ay ^emission/ {Ay) that begins when the current star for¬ 
mation rate becomes higher than the past average star 
formation rate. At lower fractions of young stars, the 
mean ratio is ^ 2.5, but this drops to ^ 2 in re¬ 
gions where the fraction of young stars is very high 
{SFRioo/{SFR) > 2). This trend is not nearly as strong 
when we limit the measurement of the current star forma¬ 
tion rate to only the past 25 Myr (not shown), indicating 
that the cumulative effect of star formation over a longer 
timescale is more important in shaping the dust emis¬ 
sion than just the short term impact of O- and massive 
B-stars. The observed trend cannot be due to contamina¬ 
tion from younger stars on the RGB discussed earlier in 
Sec. 5.3.2[ which would tend to reduce the CMD-based 
extinction and increase the plotted ratio, which is the 
opposite of what is seen in Figure [31] 

Figure also points to an interesting path fo rward. 
The recent star formation histories available in ILewisI 


et al. (2015) in principle allow the spectrum of the in¬ 


terstellar radiation field to be calculated indpendently 
of the dust emission. By assuming a vertical distribu¬ 
tion, the density of the interstellar radiation field could 
also be calculated. Including these independent measure¬ 
ments would remove a significant uncertainty in model¬ 
ing the dust emission, and would potentially lead to im¬ 
provements in modeling the dust itself. This test would 
be particularly interesting in other higher star formation 
rate intensity galaxies in the Local Group, including M33 
and the Magellanic Clouds. 


6.4. Spatially-Dependent Systematies 

The power of the CMD-based extinctions is that they 
present a direct, independent route to constraining the 
dust mass, with an angular resolution that is currently 
significantly higher than other techniques. However, 
when comparing to other ISM tracers — such as one 
may want to do to derive gas-to-dust ratios, or to con¬ 
strain the physical models of the dust — we have to be 
wary of possible systematic errors. 

We have previously discussed a number of possible sys¬ 
tematic errors that may affect the CMD-based extinc- 
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Fig. 30.— Images of the pixel-by-pixel ratio of the|Draine et al.|(|2014| emission-based Ay map divided by the resolution-matched map 
of mean extincti on (Ay) derive d from the NIR CMD, betore (lettj and after (right) applying the scale-hbias correction used in the top-right 
panel of Figure |28| (Sec. |6.2.1|). The heavy dark contours indicate where (Ay) > 1- Inspection of the uncorrected map shows that the 
mismatch between the two maps is clearly spatially dependent. The overall difference in scale in the high extinction regions falls between 
2-3, with some tendency towards higher values on the far side of the disk. The scale correction brings the agreement to roughly ±30% 
in regions with (Ay) > 1. The ratio deviates most strongly fr om 1 at low to intermediate extinctions, due to the expected systematic 
underestimates in the CMD-based extinctions (discussed in Sec. |5.3| >; the agreement in these regions is greatly improved by the correction 
of an additive bias factor, although at the cost of worse agreement at the very lowest extinctions. The correction is not as effective as in 
the high extinction regions, most likely because the bias factor would be expected to vary radially. 


tions (Sec.|4.2.5|&|5.3.2|). The majority of these possible 
systematics are due to uncertainties that are on the or¬ 
der of a few hundredths of a magnitude in the NIR color 
distribution of the unreddened RGB. While these are ex¬ 
tremely small effects in terms of the CMD, they lead to 
uncertainties of several tenths in the visual extinction, 
given that Ay = 7.56^(F110W — F160W). 

One feature of the majority of the systematic errors, 
however, is that they should only produce either smooth 
large-scale variations, or variations that correlate with 
individual WFC3/IR chips, while leaving the morphol¬ 
ogy of small scale pixel-to-pixel maps intact. One class 
of possible large scale variations we may expect is one 
that correlates with the local stellar surface density. The 
unreddened model we adopted was constructed in bins 
of local surface density to simultaneously capture pho¬ 
tometric errors, photometric depth, and s moo th stellar 
population changes with radius (see Sec. 4.2 and Fig¬ 
ure We therefore could have biases that depend on 
logi^stars if there are specific density bins in which the 
model of the unreddened CMD is incorrect. Possible ori¬ 
gins of such a bias could be: (1) the need to interpolate 
between neighboring bins when a density range contains 
few unreddened regions; (2) failure of log^^o ^stars to cor¬ 
relate perfectly with radial changes in the stellar pop¬ 
ulation, as might be expected in the rapidly changing, 
structurally complex inner regions; or (3) loss of sensi¬ 
tivity to low extinction in the high surface density inner 
galaxy where the RGB is broadened by crowding and/or 
a large metallicity spread. 

A second possible class of large scale systematic is one 
that varies with distance from the major axis. The stel¬ 
lar disk of M31 is thick (as we show in a companion pa¬ 


per), and as such there can be a range of radii present in 
projection at any gi ven position in the sky (“projection 
mixing”; Sec. 4.2.5). This effect is non-existent along 
the major axis, where the projected radius equals the 
true radius. However, the farther one moves from the 
major axis, the larger the fractional range of radii that 
are projected onto the same location. Any radial gradi¬ 
ent in the stellar population will affect the morphology 
of the RGB, and thus the superposition of a large range 
in radii can lead to offsets and spreads in the apparent 
RGB. This effect would lead to values of Ay that are 
biased in regions that are further from the major axij^ 

We can assess these two large scale effects by looking 
at the spatially-resolved correlations between the GMD- 
based extinctions and other dust tracers. Specifically, we 
break down our maps of the mean extinction into bins 
of stellar surface density and reddening fraction. The 
latter is an excellent proxy for the spread of radii present 
at a given position, because the same spread is directly 
responsible for deviations from fred = 0.5. We assign 
regions to bins of log^^Q T^stars and fred using the models 
developed when constructing the un redd e ned C MD and 
the prior on /red, respectively (Secs. 4.2 & 4.4). 

The resulting transformation of the extinction map is 
shown in Figure restricted to regions where Ay is 
likely to well-measured {{Ay) > 0.5 mag). The left hand 
plot shows the more familar map, with the locii of con¬ 
stant log]^o ^stars shown as heavy contours, and the fred 


The amplitude of the effect may also have a secondary de¬ 
pendence on azimuth in the inner regions of the galaxy, where the 
fraction of stars that are due to M31’s bulge and/or inner spheroid 
can become significant, and would not be properly captured by a 
model of a pure inclined disk. 
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Fig. 31.— The observed pixel-by-pixel ratio of the |Draine et al~] 
(|2014|) emission-based Ay map divided by the resolution-matched 
ma p of mean extinction {Ay) derived from the NIR CMD (i.e., Fig- 
urelS pF as a func t ion of the Scalo h parameter {= SFRiqq/{SFR)) 
fromlLewis et al.|(|2015|). The ratio is systematically smaller as the 
contribution trorn hot young stars increases, shown as the running 
average plotted in blue. For SFRiqq/{SFR) > 2, almost all of the 
pixels are below the globally measured value of ~2.5 (Figure [2^, 
shown as the red horizontal line. Points are color-coded according 
to the mean CMD-based extinction. 


coordinate shown as “spokes” emanating from the center 
of the galaxy. The right hand plot shows the extinction 
distribution remapped into the log^^g '^stars vs fred coor¬ 
dinate system. 

In Figure [ 33 ] we plot the pixel-by-pixel correlations be¬ 
tween the mean CMD-based extinction {Ay) (restricted 
to where {Ay) >0.5) and the scale-corrected emission- 
based extinction Ay^Qraissionl^''^‘^' Each individual cor¬ 
relation sub-panel is calculated for pixels within one of 
the bins of logir^stars versus fred shown in the right 
panel of Figure Thus, each sub-panel is equivalent 
to the upper rightpanel of Figure but restricted to a 
subregion of the galaxy. 

The correlations in Figure [33| can be used to look for 
large-scale spatially-dependent variations in the system¬ 
atic bias and in the factor of 2.53 scale factor used to 
correct Ay ^emission to a better estimate of Ay. In each 
range of log;L0 ^ stars and fred^ the diagonal line shows 
the expectation for a perfect 1:1 correlation. Data will 
fall along this line wherever Ay^emissionl‘^-^^ is a good 
estimator of Ay and there are no systematic biases in 
either extinction measurement. If the scale factor is cor¬ 
rect, but there is a systematic bias, then the data would 
have the same slope as the 1:1 line, but would be shifted 
horizontally or vertically. 

6.4.1. Spatial Variations in {Ay) Correlations 


In general. Figure suggests that the CMD- and 
emission-based extinctions are well correlated through¬ 
out the galaxy, with no strong positi on al-dep endent bi¬ 
ases like those discussed in Sec. 14. 2. 5l fc 15.3.21 There are 
two notable exceptions. 

The first is seen at very high surface densities where 
the range of emission-based extinctions is much larger 
than the CMD-based extinctions. This difference is most 
probably due to the CMD-based extinction’s lack of sen¬ 
sitivity in the very inner disk, where the CMD is too 
broad to be sensitive to any but the highest extinctions 
(Figure 10). It may also be related to an inappropriate 
model for the unreddened CMD, due to the inner region’s 
complex stellar populations, which lead to large spatial 
variations in the fractions of M31’s disk, bulge, bar, and 
inner spheroid. 

The second place where the correlation between the 
two extinction measures breaks down is at the very low¬ 
est surface densities. These regions are the origin of the 
small tail of stars with spuriously high CMD-based ex¬ 
tinctions seen in Figuremost of these are in the region 
of bad photometry in Brick 22, or in the very outskirts of 
the disk where there are both few reddened stars and low 
extinctions, which leads to systematic errors (primarily 
due to systematic offsets in the photometry across the 
WFC3/IR chip) being more noticeable. 

To quantify the trends in Figure in Figure we 
show the ratio Ay^emission/ {Ay) at {fAy) = 1, as derived 
from fitting a linear function to the correlation in each 
subpanel of Figure (i.e., as a function of log^^g '^stars 
and fr ed)^ after removing t he factor of ^2.5 correction 
to the Draine et al.| ([2014) extinctions. The plot has 
been scaled such that the overall ratio is in the middle of 
the color bar, and the range goes from zero to twice the 
observed ratio. 

Figure confirms our impressions from Figure |33l 
namely that the ratio of Ay^s is relatively uniform, ex¬ 
cept for the previously noted problems at very low sur¬ 
face densities. 


6.5. Testing the Plank Collaboration et al. ( 2 OI 4 ) 
Exetinetion Estimator 


The Plank Collaboration et al. (2014) studies of red¬ 
dening tbwardiTMnkyTVaymbleGula^^ and towards 
background QSOs strongly suggest that some of the spa¬ 
tial variation in Ay ^emission/{Ay) that we see in Fig¬ 
ures [^&[^ may be correlated with a parameter related 
to thFba ckground radiation field. In the physical dust 
models of Draine & Li (2007), the dust is heated by a dis¬ 
tribution of radiation held intensities, described by a 5- 
function at the minimum intensity Umin and a power-law 
distribution extending from Um.ir to Umax with slope a 


Plank Collaboration et al. (2014) argue that the depen- 
dence of the observed extinction on Umin suggests that 
in practice Umin is not solely tracing the minimum radi¬ 
ation field intensity, but instead is also compensating for 
variations in the properties of the underlying dust, some 
of which may be manifesting in departures from the spec¬ 
tral shape of the dust opacity expected in the Planek 857, 
545, a nd 353 bands. Altho ugh these band are not used 
in the Draine et al. (2014) fits, the SPIRE 350/im and 
500/im bands cover comparable wavelengths as the first 
two of these Planek bandpasses and both papers produce 
comparable dust mass estimates in M31 (Appendix A of 
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Fig. 32. — Maps of the mean extinction {Ay) in terms of position on the sky (left), and transformed into coordinates tracking the local 
stellar surface density and the fraction of reddened stars (right); the latter quantity is a proxy for the range of radii that fall in a given 
line of sight. Only pixels with high extinction are plotted ({Ay) > 0.5). The heavy solid lines on the left panel show lines of constant 
log^^Q Tistars, and the light diagonal lines indicate constant /red for fho same coordinate grid used in the right panel. 

no longer be robust to uncertainties in background sub¬ 
traction. The resulting fit prov ides an alternate w ay to 
estimate the extinction from the Draine & Li (2007) dust 
models: 


Plank Collaboration et al.|2014|). Thus, the|Draine et al. 
(|2014|) tits are likely to be subject to the same ettects. 


' Plank Collaboration et al. (2014) argue that Umin can 
be used to make more accurat e estimates of the e xpected 


extinction using the results of Draine & Li (2007) fitting. 
They propose an empirical linear formula to predict the 
extinction: 


^V,emission,Planck — {^Umin T X emission (1^) 

where a = 0.31 and h = 0.35 for the Milky Way’s dif¬ 
fuse gas and a = 0.33 and b = 0.27 for Milky Way 
molecular clouds. At a typical value of Umin = 0-5, 
^V,emission/^V,emission,Planck WOuld haVC ValuCS of 2.0 

and 2.3 for the diffuse ISM and molecular clouds, respec¬ 
tively. These values are comparable to the overall factor 
of 2.53 overestimate we derive here. 

We can look for a similar linear dependence of 
emission ou Umin withiu our data for M31. In Fig¬ 
ure we plot Umin versus the per c entan ge of the dust 


mass in PAHs from Draine et al. (2014), Qpah- Each 
pixel overlapping PHAl' is color-coded by its measured 
value of Ay^emission/{^v)• Wc havc cxcludcd regions 
where Ay,emission < (3-5 mag, to eliminate piexls where 
the emission and CMD-based extinctions will have the 

largest uncertainties^ _ 

Consistent with [Plank Collaboration et al. (2014), 
the 


35 


shows a tendency tor 

u 


IS 


left panel of Figure _ 

Ay^emission/(Ay) to be Overestimated where 
low. However, the value of Ay ^emission/{Ay) also ap¬ 
pears to depend on q^ah well. The right panel of Fig¬ 
ure shows a 2nd order polynomial fit to the data with 
Qpah < b%; we exclude the highest values of Qpah be¬ 
cause they are all found in a region that is very close to 
the limit where Draine et al. (2014) consider their fits to 


Ay ^emission / {^v) = 3.09 - {).2Aqpah + 0.05g|^^ 
-^AlUmin - 1 . 261 /^,, 
~\~2 AbUminQpah ^'^"^UminQpah 


(17) 


-1.261/4 


-0.181/4 


where qpah is expressed as a percentage (rather than a 
fraction). We also provide a simpler first-order linear 
correction that produces an adequate fit to the data, but 
that leaves slightly larger systematic residuals: 


(18) 


Ay ^emission /{Av) — 2.57 + O.OSqpah ~ l-24Gnm 

+0 

A7UminQpah‘ 

Figure [36] shows the resu lt of applying either the Plank 
Collaboration et al. (2014|) correction or the Umin -^Qpah 
correction trom equ.^TTp Because our map of M31 in¬ 
cludes both diffuse and molecular cl oud components, for 


the Plank Collaboration et al. (2014) correction we adopt 
a slope and intercept midway between the two very simi¬ 
lar linear corrections recommended for these two regimes 
(a = 0.32, 6 = 0.31). 

The upper panels of Figures [36] show the observed ra¬ 
tio Ay ^emission/{Ay) as a functiou of Umin^ color-codcd 
by either qpah (upper left) or the mean extinction {{Ay); 
upper right). The cyan stars show the median ratio in 
small bins of Umin^ and the magnenta line shows th e 


estimator derived in Plank Collaboration et al. (2014). 
The ho rizonta l line is the overall 2.53 scaling we derived 
in Sec. 6 .2. Ij (i.e.. Figure 28). The lower panels show 
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Fig. 33.— Correlations between the CMD-ba sed extinction and resolution-matched predictions for Ay derived from emission (after 
correct ing for Av^emission/ i^v) = 2.53; see Sec. |6.2| ), calculated in bins of log^^Q S star .s and fred: using the coordinate divisions shown in 
Figure Each subpanel contains a plot equivciTeht to the upper right of Figure with the CMD-based extinction on the horizontal 
axis arid the emission-based extinction on the vertical axis. In each subpanel, Ay is plotted between 0 and 4 mag, and data with 
^V,emission < 0.5 mag is not shown. Perfect correlations would fall along the diagonal lines. In every pixel, the correlations are strong. 
However, the re ar e subtle differences in the slopes, depending on the specific values of log^^Q stars and /red- These shifts are discussed in 

detail in Sec. 1131 
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Fig. 34.— Ratio between resolution-matched predictions for A\/ 
derived from emission (before correcting see Sec. |6.2| ) 

and the CMD-based extinction, calculated in bins of log i n '^ sto,rs 
and fred: usiug the coordinate divisions shown in Figure [3^ The 
ratio i s inf erred at (Ay) = 1 from fitting the correlations shown in 
Figure!^ but without applying the factor of 2.53 correction to the 
dust-emission based extinctions; only bins whre the data extends 
to at least (Ay) = 0.75 are shown. The color bar spans from zero 
to twice the expected ratio. The spatial variations in the ratio are 
modest for the dust-emission estimates. 

the ratio that results after correcting the emission-based 


extinction by the 2 nd-order polynomial fit from eo n. 17 
(left panel) or the Plank Collaboration et al. (2014) cor¬ 
rection (right panS)^ _ 

Figu re [36| shows that the [Plank Collaboration et alT] 
(|2014|) extinction estimator overall does a g ood job oi 
corr ecting the extinction predicted by the Draine fc] 
(2007) models (lower right). The overall scale of 
correction agrees very well, such that the ratio 
^v,emission/i^v) has a median within 25% of 1 af¬ 
ter applying the Umin correction to Ay^emission- How¬ 
ever, when compared to the result of applying the 
Umdn. + Qr^ah correct ion (lower left), the Plank Collab- 


Li 

We 


oration et al. (2014) correction leaves larger, systematic 


residuals. Both corrections leave significant scatter, how¬ 
ever, which may be due to some of the effects discussed 
in Sec. 16.31 

In Figure 1^ we compare the s patial distibution of 
Ay^p^rri .i ssion / Wv ) after making the Plank Collaboration 
et al. (2014|) correction (left) or the Umin + Qpah correc¬ 
tion from eqn. 17 (right). Th e left panel shows that the 
system atic residuals that the [Plank Collaboration et al. 
([2014[) correction leaves with Umin produce large-scale 


spatially-correl ated residuals in the predicted ex tinction 
map. With the Plank Collaboration et al. (2014) correc¬ 
tions, the extinction in the inner galaxy and half of the 
major star forming ring is overestimated, and the extinc¬ 
tion in the outer g alax y is underestimated. When corn- 


correction appears to produce larger, more spatially- 
correlated residuals than applying a single ^2.5 correc¬ 
tion factor. 

In contrast, the Umin + Qpah correction from eqn. IT7 
produces a map that appears to be free from the large 
scale systematics that were visible in Figures 
The only noticeable exceptions are the two re gion s of 
high star formation intensity discussed in Sec. (6.3[ (see 
Figure 31) where the corrected value of Ay^^mission is 
too low, and a region near the end of M3Fs long bar in 
the inner disk, where the corrected value of Ay^^mission 
is too high. This latter mismatch may be due to an 
error in {Ay) associated with rapidly changing stellar 
poplations at this position. Other than these regions, 
and the expected convergence of the uncorrected value 
of Ay ^emission and {Ay) in regions where Ay ^ 0 in the 
outermost disk, the corrected map looks quite uni form. 
This suggests that our earlier suspicion (Sec. 6.3) that 
some of the observed systematic variations were tracing 
the variation in Qpah was correct. 

6.5.1. Interpreting the Dependenee on qpah 


As Plank Collaboration et al. (2014) noted, the depen- 
dence of Ay pm.issi nni{Ay) on othcr parameters of the 


pared to Figure 30, the Plank Collaboration et al. (2014) 


Draine & Li (2007) fits does not necessarily suggest that 
the true values of those parameters are actually chang¬ 
ing. Any fitting procedure will tend to produce corre¬ 
lated residuals in the values of parameters when the un¬ 
derlying model is not perfect. The mass contained in 
PAHs is quite small (typically less than 5%), and they 
are responsible for very little of the extinction in the 
NIR. Thus, the fact that we see a correlation between 
the Ay ^emission! {Ay) and qpah is unlikely to indicate that 
variations in the PAH fraction are actually causing vari¬ 
ations in the extinction. Instead, it is likely that both 
Umin and qpah are compensating for subtle differences 
in the spectral shape of the dust’s opacity and emissiv- 
ity, and for complexities in how the interstellar radiation 
field relates to dust emission. 

One interesting feature in Figure points to possi¬ 
ble connections to the role of dust heating and to the 
differences between dust in the diffuse ISM and in star 
forming regions. The upper panels of Figure ap¬ 
pear to show two broad diagonal sequences that have 
a similar variation of Ay ^emission/{Ay) with Umim but 
that are shifted horizontally from each other. At low 
values of Umin < 0-2, the majority of points have low 
PAH fractions and low extinctions. At higher values of 
Umin >0.2, the sequence is dominated by regions with 
higher PAH fractions and higher extinctions. This bifur¬ 
cation may reflect that the appropriate conversion be¬ 
tween emission and extinction largely separates into two 
different regimes - one associated with the diffuse ISM 
and one associated with higher star formation rate in- 
tensities and large r typical gas columns. Although the 
Draine & Li (2007) models attempt to include this effect 
by allowing some fraction of the dust to be heated by 
harder radiation fields associated with photo-dissociation 
regions (PDRs), it is unlikely that a single physical dust 
model can fully capture the complexity of dust heating 
and composition that is found in a typical PDR. 

In light of this, although Figures [35 H 37 | show that the 
correction to Ay ^emission given in eqn. f^does an excel¬ 
lent job of predicting the extinction in JH31, the correc- 
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Fig. 35. — The ratio between the emission-based extinction and the mean CMD-based extinction, as a function of Umin and qpah derived 
in|Draine et al.| (|2014|). The left panel shows the data and the right panel shows the results of a 2nd order polynomial fit. Th e data show 
that tJie observed ratio depends on both Umin and Qpah^ not just on Umin as was assumed in [Plank Collaboration et al.|(| 2014]). The range 
of the color scale is chosen such that the grey value indicates the mean ratio of Av^emission/{^v) = ^-53 found in Sec^6.2.1| The data 
is only plotted for points with Ay^emission > 0-5, where uncertainties and biases in both methods are sm aller. The fit was re stricted to 
regions where Qpah < 6%, because the regions with higher qpah are all found at the outer limits of where [Draine et al.j (|2014| considers 
their fits to be reliable. For clarity, small random offsets have been added to both Umin and qpah fo compensate for the quantization of 
the original data and to better show the individual values of Ay^emission/{Ay)- 


tion should not necessarily be expected to do an equally 
good job in other circumstances. The joint distribution 
of Umin and Qpah IS likely to reflect the dust heating and 
composition that is particular to M31, to first order. The 
second order dependence of Ay^emission/i^v) on Umin 
and Qpah seems likely to arise front devia tions between 
the assumptions of the|Draine & Li|(|2007|) physical dust 
model and the real benavior of IVldl's dusty ISM. Nei¬ 
ther of these situations will necessarily be reproduced in 
galaxies with very different star formation rates, metal- 
licities, or dense gas fractions. Likewise, if the variation 
of Ay^^Qmission I {Ay) with Umin and Qpah Is duC tO the 
fit compensating for subtle mismatches in the underly¬ 
ing dust model, then fits that use data taken in differ¬ 
ent bandpasses and/or with different physical resolutions 
may not compensate for these effects in exactly the same 
way. As a result, we do not expect that the correction in 
eqn. 17 will be one that holds generically in all situations. 
The overall scaling of a factor of ^2 — 2.5 seems robust, 
given that comparable corrections were needed both in 
M31 and in the diffuse and the molecular components 


of the Milky Way (jPlank Collaboration et al.||20i4| 
the exact second-order dependence of Ay ^emission/ 
on Umin and Qpah may be variable. 


but 

(Av) 


7. CONCLUSIONS 

In this paper we have developed a new method for map¬ 
ping extinction in external galaxies. The method is based 
on using NIR observations of RGB stars. These stars 
have a very narrow distribution in NIR color as a func¬ 
tion of magnitude, allowing one to diagnose extinction as 
a redward shift in the position of the RGB on a GMD. We 
model the NIR GMD as a combination of a foreground 
of unreddened stars, and a background of stars that he 


behind a thin layer of dust. The background stars sam¬ 
ple the distribution of extinctions within the dust layer, 
which we characterize with a single log-normal function. 
We fit for the fraction of stars behind the dust layer, and 
the characteristic parameters of the log-normal (median 
extinction, dimensionless width), from which we derive 
the mean and dispersion of the distribution of Ay. We 
carry out this analysis using stars from the PHAT sur¬ 
vey of M31, selected in 25 pc (6.6") bins with 12.5 pc 
sampling. 

The resulting map is the highest resolution tracer of ex¬ 
tinction available in M31 to date. More importantly, it 
offers a completely independent assessment of the prop¬ 
erties of the cold ISM, allowing it to be compared with 
other widely-used tracers, including emission from dust 
and from cool gas (Hi and GO). 

Our maps reveal a factor of ^2.5 offset between the 
extinction directly measured from RGB stars and the 
extinction inferred from rece nt maps of the dust emis¬ 
sion (e.g., Draine et al.||2014|). This difference is unlikely 
to be due to the INIR (JMD fitting, and instead suggests 
that there is significantly more emission associated with 
a given amount of dust extinc tion than expected f or the 
most recent calibration of the Draine & Li (2007) mod¬ 
els. Given how important these models have become to 
the study of external galaxies, understanding the origin 
of this offset is vital. Similar offsets have been seen in 
the Milky Way (jPlank Collaboration et al.||2014|), where 
the y were found to hav e a second order depenaence on 
the Draine & Li ([2007) estimates of the strength of the 
background radiation field. We show that the offset in 
M31 depends on the measured percentage of PAH’s as 
well. We provide a fitting formula to allow the extinc¬ 
tion to be estimated from the results of fitting to the 
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Fig. 36.— The observed ratio between the emission-based extinction and the mean CMD-based extinction, as a function of Umini color 
coded by Qpah (fop left, bottom row) or {Ay) (top right). Large cyan s tars indicate the median value of A y^emissionl(Ay) in bins of Umin- 
In the top row, the solid magnenta line is the ex pectation based on the|Plank Collaboration et al.|(|2014|) correction, and the horizontal red 
line is the mean correction derived in Sec. 16. 2 . ^ The bottom row shows tJie ratio atter correcting ^y^ernission using the joint Umin + Qpah 
correction shown in Figure [S^ (left) and tne|Plank Collaboration et al.|(|2014|) correction based on Umin (right). The correction based on 
both Umin and Qpah shows much weaker systematic residuals, except at tJie very lowest values of Umin- The spatial correlation of the 
residuals is also reduced (see Figure [3^ . For clarity, small random offsets have been added to Umin to compensate for the quantization of 
the original data and to better show rne individual values of Ay^^mission/(Ay). 


Draine & Li (2007) models. 

Given the success of this technique, it is natural to 
think of possible extensions. Within M31, one might 
consider including optical data in the analysis to provide 
better sensitivity at low extinction. However, unlike the 
NIR, the color of the optical RGB is quite sensitive to 
age, leading to a much broader RGB for galaxies with 
extended star formation histories, and greater sensitivity 
to stellar population gradients. Given that the efficacy 
of the technique depends on having sufficient reddening 
to move stars out of the unreddened GMD, it is not im¬ 


mediately clear that optical data’s greater sensitivity to 
reddening would counteract the problem of an intrinsi¬ 
cally broader RGB. We address these issues more quan¬ 
titatively in an appendix (Sec. 0 - 

Some of the difficulty in incorporating optical and 
UV data into an RGB-only analysis can be sidestepped 
by carrying out joint UV+optical+NIR fitting of every 
star’s spectral energy distribution (Gordon et al 2015, 


the IVlagef 


anic Glouds, but used UV and optical data 


Zarit- 


submitted). Thi s ap proach was very succe ssful in 
sky et al.| (2002) and Zaritsky et ay(2004|)’s analyses of 
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Fig. 37. — Maps of the ratio between the emission-based extinction and the mean CMD-based extinction, after correcting Av^emssion by 
applying either thelPla nk Collaboration et al.|(|2014|> co rrection that uses only Umin (left) or a polynomial correction based on both Umin 
and Qpah (right), ine ([Flank (Jollaboration et al.|2Dl4|) correction does a good job of correcting the overall normalization, but introduces 
a clear radial bias and does not correct some of tne largest, spatially-coherent features in the map. The correction based on both Umin 
and qpah produces a much more uniform map. The notable discrepancies are in t he i nner disk where (Ay) values are known to be more 
uncertain, and in the two most actively starforming complexes identified in Figure [sT] 


alone, as was appropriate for the much lower extinctions 
in these galaxies. Fitting PHAT’s 6-filter photometry 
should offer an excellent opportunity to probe down to 
lower extinctions and to assess the detailed connections 
between dust and different ages of stellar populations, 
although with more dependence on the accuracy of the 
underlying stellar models. 

As another possible extension, one might also think of 
applying this technique to other galaxies, with M33 being 
an obvious target. In more distant galaxies, however, the 
apparent width of the RGB will be a limiting factor when 
using current facilities; more distant galaxies have higher 
levels of crowding that limit the photomet ric accuracy 


possi ble with a given telescope+camera (e.g. Olsen et al. 


2003), leading to a wider RGB. However, a next gen¬ 
eration telescope (either 30m class from the ground, or 
an 10-12m in space) with a high resolution NIR camera 
would make this technique viable throughout the local 
volume. 
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APPENDIX 

APPLYING THIS TECHNIQUE AT OTHER WAVELENGTHS 

Although we have developed and applied this technique in the NIR, it is natural to consider whether it might also 
be useful in other wavelength regimes. In particular, by incorporating optical data we could potentially improve our 
sensitivity to low extinctions, because of the greater impact of reddening and extinction at bluer wavelengths. As we 
now show, the increased sensitivity to reddening comes at a cost, because the underlying CMD is far more sensitive 
to population effects in bluer filter. 


,ian/Av) after U^i„+qpah correction to 
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The key issues that affect the sensitivity of our technique are: (1) the amount of reddening produced by a given 
Ay; (2) the amount of reddening needed to move a star clearly out of the RGB; (3) the spatial uniformity of the 
unreddened RGB; (4) the maximum Ay that can be detected before losing a star below the photometric limit; and 
(5) the total number of well-measured stars on the RGB. Incorporating data from bluer wavelengths improves the first 
issue, but its effect on the remaining issues is unclear. 

To assess the impact of these issues outside the NIR, we construct GMDs of the unreddened RGB in a variety of 
filter combinations. In Figure 38 we show Hess diagra ms constructed with in Brick 12, Field 7, which has extremely low 
reddening as judged from boththe NIR GMD and the Draine et al. (2014) maps. To allow us to compare arbitrary filter 
combinations, we use the simultane ous six filter * . gst phot ometry described in Williams et al (2014), which supersedes 

2012); these photometry have not had cuts applied to optimize 
GB, but are sufficient for our purposes here. The horizontal 


W 


the first generation photometry in Dalcanton et al. 
rejection of spurious measurements redward of the 
purple bars indicate the median photometric uncertainty at 5 different magnitudes; the uncertainties are a quadrature 
sum of the photon coun ting uncertainties reported by DOLPHOT and a fiat 0.015 mag uncertainty due to calibration 
(Dalcanton et al. 2012). The diagonal red line shows the effect of AAy = 2mag, using Ap^j^w= 1.1591 and 
Af 814 W/Ay = O.o962, which we derived analogously to Apw^y/ / Ay and ^fi6QWMv above in Sec. 

The first, and most important, figure of merit needed to interpret Figure is the relative wid^th of the RGB 
compared to the length of the reddening vector. If the RGB is narrow compared to the reddening vector, then even 
modest Ay will move a star clearly off the unreddened sequence, increasing the sensitivity to low dust columns. We 
facilitate this comparison by plotting each GMD with a color range that corresponds to a constant AAy, such that 
the reddening vector has the same width in each panel. With this display choice, GMDs that appear wider are indeed 
“wider compared to a given AHy”. 

Inspection of Figure [38| shows that the RGB is intrinsically narrower in the NIR than in any other filter combination. 
In other words, althou^ the reddening in FllAW — F1601F is smallest for a given Ay, our ability to actually 
detect that reddening is largest. In terms of Ay, the NIR RGB has a width of AAy = 0.38 mag for stars with 
19.5 < F1601F < 20.5, whereas the optical GMD (upper left) has AAy = 0.54 mag. 

The narrowness of the NIR GMD comes in large part because of the cancellation of age and metallicity in RGB 
isochrones. Older stars tends to have redder RGBs, but these stars also tend to have lower metallicities, which drives 
their colors back blueward. In the NIR, these effects are intrinsically small, and also largely cancel out. In contrast, in 
the optical these effects are much larger and do not effectively cancel out. This unfortunate fact leads to broader GMDs 
in the optical. Moreover, it points to increased sensitivity of the optical GMD to variations in the underlying stellar 
population’s age and/or metallicity. While this is good for stellar population studies, it means that one expects much 
larger spatial gradients in the properties of the unreddened RGB, making the construction of an accurate, smoothly 
varying “unreddened” GMD far more difficult. Thus, by the criteria listed above, issues (2) and (3) are negatively 
affected by including bluer filters. 

Issue (4) is also further compromised when moving to bluer filters. Inspection of the optical GMDs in the top row 
shows that the blue magnitude limit in F475IF cuts off stars with more than Ay ^ 2. In contrast, the redder GMDs 
in the bottom row allow a much larger range of Ay to be detected, giving better constraints on the median Ay and 
the shape of the reddening distribution. This limitation can not necessarily be overcome by a larger investment of 
telescope time, given that the data is already crowding limited. 

The one area where switching to optical filter may improve sensitivity is in the number of stars (i.e., issue #5). 
Comparing photometric errors (magenta bars) to the width of the RGB shows that bluer filters can accurately track 
the RGB down to fainter GMD features, increasing the number of stars that can be used in a given area. This difference 
potentially allows one to increase the resolution of the extinction maps. However, the higher stellar density would have 
only modest effects in high extinction regions, since most of the newly detected stars are close to the magnitude limit. 

Rather than switching entirely to the optical, it may be that some gains could be achieved by using F814IT —F160IF 
to measure the reddening. The RGB is less than 5% wider than in the NIR (in terms of AAy), so the sensitivity to 
small extinctions should be comparable. However, the photometric errors are proportionally smaller, due to the wider 
color baseline and the better accuracy in FS14W compared to FllOW. These improved errors should allow us to use 
stars further down the RGB, which increases the spatial density of stars. One can then use smaller analysis pixels to 
increase the spatial resolution of the extinction map. Unlike in the optical, these additional stars are comfortably above 
the magnitude limit, even for significant reddening, and thus they should provide useful leverage on the distribution 
of reddenings. It is not clear, however, whether these gains trump the likely increase in the susceptibility of the RGB 
structure to age and metallicity variations. 


SMOOTH MODELS FOR THE DUST DISTRIBUTION 

During the course of this work, we also considered a model where the reddening seen in the GMD was produced by 
a smooth layer of dust embedded within the stellar disk. As we show here, this assumption produces distributions of 
extinction that are inconsistent with the data, leading to the need for the alternative “thin screen with a log-normal 
distribution of extinction” model adopted in this paper. 

We calculate the expected distribution of stellar extinctions for a smooth dust model as follows. We assume that 
the dust is embedded within a stellar disk with a scale height We model the stellar disk using the probability 
Pi,{z; z^) that a star is found at height 2 : relative to the midplane of the stellar disk (assuming that positive values of z 
correspond to stars that are closer to us than the midplane at 2 : = 0). The stellar density distribution can be assumed 
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Fig. 38. — Comparing the sensitivity of different filter combinations to reddening and extinction. Plots show logarithmically-scaled Hess 
diagrams of the RGB region in 4 different filter combinations, for one of the lowest reddening regions of the PH AT survey (Brick 12, Field 
1, plotting the ‘ ‘gst’ ’ photometry). The horizontal magenta bars indicate t he median reported ph otometric uncertainty at the plotted 
magnitude, including an additional systematic component of cf color = 0.015 (see|Dalcanton et al.|2012|); these bars do not include errors due 
to crowding, which will increase significantly near the limiting magnitude of the data. In general, tlie width of the RGB is significantly larger 
than the photometric uncertainties, indicating that variations in the age and metallicity of the stellar population is primarily responsible 
for broadening the RGB. This effect is worse in the optical filters, but negligible in the NIR, where the RGB is intrinsically narrow and has 
a width almost entirely dominated by photometric uncertainties. The red bar shows the shift in color and magnitude that would be caused 
by /S^Ay = 2 magnitudes of extinction; and the range of plotted colors is adjusted so that it covers the amount of reddening caused by a 
constant AAy = 10 magnitudes of extinction. The intrinsic width of the RGB in terms of /S^Ay is much broader in the optical filters than 
the NIR, which cancels out any possible benefit which may have resulted from the optical’s greater sensitivity to extinction, compared to 
the NIR. 


to be a Gaussian, an exponential, or any other reasonable functional form. 

We first assume that the dust is smoothly distributed within a vertical distribution Pd{z; Zd^ zq) with scale height 
Zd, centered at a height zq relative to the midplane of the stellar disk. We assume that the total extinction of the 
integrated dust column is Ay^tot- The integrated extinction seen for a star at height z is then 

POO 

Av{z)=Av,tot J Pd{z\Zd,zo)dz. (Bl) 

If we assume that the vertical distribution of dust is an exponential and let z' = z — zq, then 


Av{z = z' A- Zq) 


W,tot 


2;' - 2;' 


+ 
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Fig. 39. — Probability distributions for a smooth exponential dust layer with scale height z^ust embedded within a Gaussian stellar disk 
of scale height Zstars: for different ratios between the scale heights of the dust and stars. When the dust layer is thin compared to the 
stars, the stars have a bimodal distribution of Ay, falling in an either unreddened or reddened peak. This bimodality becomes less distinct 
as the thickness of the dust and stars become comparable. The left hand panel shows a dust layer centered at the stellar midplane and the 
right hand panel shows a dust layer that is shifted to the far side of the stellar midplane by one half the scale height of the stars. When 
the dust lies on the far side of the midplane, a larger fraction of stars have low extinctions. 


The probability pa{Ay) dAy of finding a star with an extinction between Ay and Ay + dAy is then 


PA{Ay) dAy =p^{z{Ay); Z^) 
= p^{z{Ay);z^) 


dz 


dAy 


dA\ 


Zd 


Av,totl‘2 ^ 

p'l Av,totl^ b'l 


7 dAy 


where z' = z{Ay) — zq and z{Ay) is the inverse of equation |B2| 


/ A \ ^Ay 

z{Ay) = Zo- Zd In 


|AAy| - AA 
lA^vl 


Al 


AAi 


\AAy\ yAy^tot/^^ 


(B3) 

(B4) 

(B5) 


assuming A Ay = Ay^totl‘^ — 

The resulting distribution of pA{Ay) is shown in Figure 39 for the case of a Gaussian distribution of stars. When 
the scale height of the dust is less than that of the stars, tne probability distribution has two strong peaks — one 
at Ay = 0, corresponding to unreddened stars seen on the near side of the dust layer, and one at Ay = Ay^toti 
corresponding to stars seen on the far side of the dust layer, after they have experienced the entire dust column. The 
relative height of the two peaks depends on the position of the dust layer relative to the stars, with the unreddened 
peak becoming stronger when the dust layer is on the far side of the stellar midplane. As the scale height of the dust 
increases, there are fewer stars with either zero reddening or the maximum reddening (i.e., all stars are embedded 
within the dust), and the two peaks move together towards a typical extinction of half the total integrate d ex tinction. 

To transfer pA{Ay) to the observed reddening distribution for the RGB, one must convolve equation B3 with the 
intrinsic distribution of RGB colors along the reddening vector. For an intrinsic Gaussian distribution of RGB colors 
and a thin dust layer, one would then expect to see two skewed Gaussian peaks - an unreddened Gaussian skewed 
towards redder colors, and a highly reddened Gaussian skewed toward bluer colors - provided that the expected total 
reddening from dust is significantly larger than the intrinsic width of the RGB colors. The reddened and unreddened 
Gaussians should have identical widths when the dust is centered on the mid^ane, but can have somewhat different 
widths if the dust lane is not aligned with the stars (right hand panel of Figure 39), in the sense that the peak containing 
a higher fraction of stars should also be wider. However, there are many regions in our data which appear to violate 
this condition. Frequently the reddened distribution is nearly an order of magnitude broader than the unreddened 
peak, which cannot be accounted for in models with a smooth dust layer. We therefore reject a smooth dust layer as 
being an appropriate model for these observations. 
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